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ABSTRACT  OF  THESIS 


MODEL  SIMULATION  OF  A  LOCALIZED  HIGH  INTENSITY  HEAT  SOURCE 
INTERACTING  WITH  COOLED  METAL  PLATES 


The  basic,  generic  problem  of  a  localized  high  Intensity  beat 
source  directed  against  one  surface  of  a  plate  of  finite  thickness  was 
investigated  using  the  finite  element  program  ANSYS .by -Swanson  Analysis'* 
Systems ,  -Inc .  After  reviewing  similar  work  in  nuclear  fuel  and  laser 
machining,  ANSYS  was  verified  against  a  known  solution.  ANSYS  was  used 
to  create  a  model  that  yields  minimum  heat  transfer  coefficients  needed 
to  prevent  the  initiation  of  melting  in  thin  aluminum,  titanium,  and 
stainless  steel  (AISI  304)  plates.  These  heat  transfer  coefficients 
were  converted  into  minimum  local  Nusselt  numbers  and  graphed  against 
local  Nusselt  number  correlations  for  constant  temperature  flat  plates 
in  forced  and  free  convection  regimes.  A  detailed  listing  of  both 
laminar  and  turbulent  correlations  is  presented  along  with  references. 
The  suitability  of  liquid  sodium,  air,  and  water  (under  high  pressure) 
as  coolants  for  a  source  intensity  of  2 . 0  ($<  107  W/m^was  examined.  For 
free  convection,  only  liquid  sodium  cooling  a  titanium  plate  is 
feasible.  For  forced  convection,  liquid  sodium  is  feasible  in  laminar 
flow  for  all  three  plates  with  velocities  ranging  from  0.28  m/s  to  1.09 
m/s ,  Water  is  feasible  for  aluminum  and  titanium  in  turbulent  flow  at 
velocities  of  approximately  4  m/s. 
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Chapter  1 
INTRODUCTION 

The  objective  of  this  work  is  to  investigate  the  basic  problem  of  a 

localized  high  intensity  heat  source  directed  against  one  surface  of  a 

plate  of  finite  thickness  and  cooled  on  the  other  surface  by  free  or 

forced  convection.  The  applicability  of  this  problem  extends  across  a 

wide  variety  of  fields .  In  the  nuclear  industry  this  situation  is  known 

as  the  'hot-spot'  problem  and  is  encountered  in  various  fuel  assemblies, 

particularly  fuel  pellets.  In  industries  that  use  laser  machining,  the 

applications  include  possible  cooling  problems  encountered  during 

cutting,  drilling,  or  welding.  In  the  defense  industry,  with  the  advent 

of  lasers  as  viable  weapons,  the  problem  extends  to  the  development  of 

possible  defensive  measures  incorporating  various  cooling  systems.  In 

the  space  industry,  possible  applications  include  the  minimization  of 

laser  effects  on  satellites  as  well  as  the  study  of  'hot  spot' 

development  during  the  burning  of  solid  rocket  fuel.  Other  possible 

applications  include:  cooling  requirements  of  energy  collectors  for 

solar  power;  the  study  of  cooling  effects  on  electrical  arcing  and 

conductor  diameter  variation;  and  local  heat  generation  in  gun  barrels 

due  to  raised  spots.  Additionally,  the  basic  problem  in  all  these 

fields  can  include  a  localized  high  intensity  heat  source  superimposed 

over  an  already  existing  heat  flux  or  temperature  field.  In  any  case, 

all  of  the  above  examples  include  the  presence  of  a  localized  high 

intensity  heat  source  where  cooling  requirements  must  be  determined. 

Thus,  a  computer  model  that  addresses  a  part  of  this  basic  problem  would 

prove  very  useful.  Such  a  model  is  presented  here,  with  the  above 

applications  simplified  to  examining  the  case  of  a  thin  vertical  metal 

1 


2 


I 

plate  subjected  to  steady- state  heating  conditions.  i 

The  specific  question  this  computer  model  is  designed  to  answer  is: 
what  minimum  heat  transfer  coefficient  is  needed  to  prevent  a  localized  \ 

f 

high  intensity  heat  source  from  initiating  melting  on  the  front  face  of 
a  metal  plate  assuming  only  the  back  face  of  the  plate  is  cooled?  The 
data  assumed  known  is:  the  intensity  of  the  heat  source;  the  area  of 
the  front  face  of  the  plate  affected  by  the  heat  source;  the  thickness 
of  the  plate;  and  the  thermal  properties  of  the  plate.  A  proposed  heat 
transfer  coefficient  is  supplied  by  the  user,  and  the  model  then 
provides  the  top  ten  temperatures  (based  on  node  location  from  a  finite 
element  program)  on  the  front  and  back  faces  of  the  plate  in  order  to 
determine  if  melting  would  occur  and  to  analyze  cooling  conditions .  The 
user  continues  to  adjust  the  value  of  the  proposed  heat  transfer 
coefficient  until  the  maximum  temperature  on  the  front  face  is  just 
below  the  melting  point  of  the  material  under  examination.  Since  the 
front  surface  of  the  plate  outside  of  the  affected  area  of  the  heat 
source  is  considered  to  be  adiabatic,  Including  the  plate  edges,  the 
resulting  heat  transfer  coefficient  represents  an  upper-bound  value  of 
the  minimum  cooling  required  to  prevent  the  initiation  of  melting  on  the 
front  face.  This  resulting  minimum  heat  transfer  coefficient  is  used  to 
determine  the  feasibility  of  potential  coolants.  The  latitude  available 
to  the  user  in  varying  the  material  of  the  plate  is  limited  only  by  his 
knowledge  of  its  thermal  properties.  The  limitations  on  the  dimensions 
of  the  plate  were  not  fully  investigated.  A  0.05  m  square  plate  was 
used  in  the  testing  of  the  model  with  a  thickness  of  0.0004  m.  The 
model  was  verified  for  thicknesses  of  up  to  0.0032  m  in  the  case  of 
aluminum.  Further  investigation  was  precluded  by  time  constraints. 
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The  two  main  fields  that  generated  interest  in  and  contributed  most 
to  the  development  of  this  model  are  nuclear  fuel  and  laser  machining. 
Thus  their  applications  of  the  model  will  be  explored  in  detail.  In 
nuclear  fuel,  the  effect  of  ’hot  spots'  in  fuel  elements  can  be  modelled 
by  estimating  worst-case  heat  intensities  from  contact  of  the  fuel  with 
the  cladding.  (This  contact  results  from  pellet  cracking  and  relocation 
as  well  as  fuel  deformation  that  occur  during  normal  operation  of  the 
reactor.)  These  intensity  estimates  plus  representative  values  of 
cooling- fluid  flow  and  cladding  thickness  can  then  be  entered  into  the 
model.  To  provide  this  information  to  the  program  requires  experimental 
or  calculated  data  of  typical  heat  intensities  due  to  ’hot  spots'  in 
fuel  elements.  Fortunately,  much  data  applicable  to  this  problem  has 
already  been  published  in  the  form  of  contact  conductance  values  [1-7]. 
Worst-case  heat  intensities  can  then  be  calculated  by  multiplying  these 
experimental  contact  conductances  by  estimated  worst-case  temperature 
distributions  obtained  from  either  accident  condition  studies  or  normal 
operation  studies.  Characteristic  temperatures  can  thus  be  either 
calculated  from  current  reactor  models  [8-12]  or  obtained  directly  from 
test  data  [12-14].  Thus,  with  these  data,  the  computer  model  can 
determine  the  new  cooling  conditions  necessary  to  prevent  the  cladding 
from  melting  due  to  ’hot  spots'  in  the  pellets. 

In  laser  machining,  the  model  is  especially  applicable  to  laser 
drilling.  By  assuming  free-convection  cooling  on  the  backside  of  the 
plate,  the  model  takes  on  the  form  of  a  laser- drilling  problem  with  the 
localized  high  intensity  heat  source  representing  the  laser  beam.  The 
model  does  not  take  into  account  the  effect  of  laser  induced  plasma;  but 
as  noted  by  Herziger  [IS],  the  undesirable  properties  of  the  laser 


induced  plasma  cloud  can  be  avoided  (or  at  least  minimized)  by  using  an 
appropriate  pulse  program  for  pulsed  lasers.  For  continuous  wave  (cw) 
lasers,  experiments  demonstrate  that  machining  with  them  is  only  of 
advantage  in  a  few  cases,  preferably  when  the  generation  of  laser 
induced  plasma  can  be  avoided.  So  neglecting  the  effects  of  laser 
induced  plasma  on  the  model  does  not  invalidate  its  usefulness  to  the 
area  of  laser  drilling  and  at  the  same  time  maintains  the  model's 
generic  nature. 

In  the  following  study,  a  review  of  some  of  the  similar  work  done 

in  nuclear  fuel  and  laser  machining  is  presented.  Then  the  rationale 

is  given  for  choosing  the  finite  element  package,  ANSYS,  by  Swanson 

Analysis  Systems,  Inc.  as  the  basis  for  the  model.  The  basic  principles 

behind  the  operating  procedures  of  ANSYS  are  summarized  and  the  program 

is  verified  by  a  comparison  run  against  analytical  data  provided  by 

Torvik  [16].  Though  the  basic  problem  being  investigated  includes  flat 

plates,  walls  of  cylinders,  and  walls  of  spheres,  only  the  case  of  the 

vertical  flat  plate  is  examined.  Aluminum,  titanium,  and  AISI  304 

stainless  steel  plates  are  examined  under  steady-state  heating 
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conditions  in  the  1  x  10  to  6  x  10  W/m  range  in  order  to  determine 
the  upper  bounds  of  the  minimum  cooling  requirements.  The  heat  transfer 
coefficients  obtained  from  the  model  are  used,  along  with  local  Nusselt 
number  correlations  for  free  and  forced  convection,  to  determine  the 
feasibility  of  air,  water,  or  sodium  meeting  the  minimum  cooling 
requirements  for  the  materials  and  conditions  examined.  Relationships 
of  heat  transfer  coefficient  as  a  function  of  heat  source  intensity, 
area  of  plate  affected  by  the  heat  source,  and  conductivity  are 
presented  in  graphical  and  tabular  form.  Heat  transfer  coefficient  as  a 


Chapter  2 
BACKGROUND 


2.1  Introduction 

The  generic  problem  of  a  localized  high  intensity  heat  source 

directed  against  a  vertical  flat  plate  with  either  free  or  forced 

convection  cooling  on  the  opposite  side  has  received  only  partial 

coverage  in  the  literature.  In  fact,  two  separate  literature  searches-- 

one  of  nuclear  fuel  and  one  of  laser  machining- -produced  no  direct 

handling  of  the  problem.  The  literature  search  of  nuclear  fuel  involved 

the  DOE  Energy  Database  extending  back  through  1974.  The  search  of 

laser  machining  involved  the  DOE  Energy  Database  extending  back  through 

1974;  the  Metadex  Database  back  through  1966;  the  Compendex  Database 

back  through  1970;  and  the  Dissertation  Abstracts  Online  Database  back 

through  1961.  The  work  done  in  the  laser  material  processing  field 

deals  almost  exclusively  with  a  localized  high  intensity  heat  source 

directed  on  a  flat  plate  (vertical  or  horizontal),  but  cooling  on  the 

opposite  side  of  the  plate  is  not  considered  (backside  is  either 

infinite  or  adiabatic).  In  the  nuclear  field,  most  of  the  work  centers 

around  determining  the  conductance  of  the  gap  between  the  fuel  pellet 

and  its  cladding.  'Hot  spots'  arising  from  relocation  of  the  fuel  in  the 

pellet  due  to  normal  operation  or  fuel  deformation  are  treated  in  either 

a  statistical  manner  or  through  the  determination  of  a  contact 

conductance.  Clad-to-coolant  heat  transfer  is  generally  calculated  by 

use  of  correlations  for  film  coefficients  in  subroutines  of  the  various 

computer  models  of  reactor  operation.  Thus  cooling  is  considered  in  a 

generalized,  lumped  correlation.  Melting  is  considered  only  as  it 

applies  to  the  fuel  since  the  fuel  is  generally  at  much  higher 
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temperatures  than  the  cooled  cladding.  A  summary  of  representative  work 
in  both  fields  follows. 


2.2  Representative  Work  in  Nuclear  Fuel 

Gamier  and  Begej  [17,  18]  come  close  to  the  generic  problem  in 
their  study  of  thermal  gap  and  contact  conductance  between  depleted 
uranium  dioxide  (UC^)  and  Zircaloy-4  (Zr4)  using  the  Modified  Pulse 
Design  (MPD)  technique.  MPD  is  a  transient  technique  employing  a  laser 
heat  pulse  and  a  signal  detector  to  monitor  the  thermal  energy 
transmitted  through  a  UOg/Zr^  sample  pair  which  are  either  in  contact  or 
separated.  This  technique  is  a  modification  of  the  classical  heat  pulse 
(flash)  method  [19]  commonly  used  to  determine  the  thermal  diffusivity 
and  conductivity  of  materials.  MPD  uses  a  laser  to  supply  a  heat  pulse 
to  the  front  face  of  a  IK^-Zr^  sample  pair  that  are  in  light  contact  or 
slightly  separated.  An  intrinsic  fast  response  thermocouple  monitors 
the  subsequent  rise  on  the  back  face  of  the  Zircaloy.  Using  this 
transient  temperature  response  and  knowledge  of  the  thermal  properties 
of  the  specimens,  the  thermal  gap  conductance  is  determined  as  a 
function  of  changing  variables  (i.e.  temperature,  mean-plane  gap  width, 
gas  composition,  etc.).  The  mathematical  description  of  MPD  is  very 
similar  to  that  of  the  present  work;  however,  several  of  its  basic 
assumptions  are  quite  restrictive.  In  MPD,  only  one -dimensional  heat 
conduction  is  considered,  due  to  the  experimental  assembly  construction. 
Thermal  properties  of  the  specimens  are  considered  constant  since  the 
temperature  transients  are  no  greater  than  4°K.  The  unknown  quantities 
are  the  heat  losses  of  the  front  and  back  faces  and  the  heat  transfer 
across  the  gap.  Laplace  Transforms  are  used  on  the  governing  equations 
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and  the  very  short  laser  pulse  is  considered  to  be  a  Dirac  Delta 
function.  A  method  developed  by  Peckham  [20],  that  attempts  to  minimize 
the  s’  i i  of  the  squared  errors,  is  chosen  for  determining  the  unknown 
parameters.  Gamier  and  Begej  conclude  that  the  MPD  technique  is 
reliable  for  gap  conductance  measurements  but  not  for  solid-solid 
contact  conductance  due  to  the  difficulty  in  applying  suitable  contact 
pressures  to  the  specimens.  They  also  provide  a  Fortran  program  for 
calculating  the  best  fitting  convective  coefficients  for  a  given  set  of 
data. 

Suzuki  et.  al.  [21]  present  a  code  which  calculates  a  cooling  flow 
distribution  in  a  reactor  vessel  and  couples  it  with  a  code  which 
analyzes  temperature  distribution  in  a  fuel  assembly.  Various 
analytical  models  for  coolant  flow  distribution  in  the  reactor  vessel, 
for  pressure  drop  calculations,  and  for  temperature  distribution 
calculations  in  coolant  and  fuel  elements  in  the  assembly  are  presented 
in  detail.  An  expression  is  also  given  for  calculating  'hot  spot' 
temperatures  but  only  through  the  use  of  statistical  uncertainties  in 
nominal  temperature  differences.  This  allows  the  'hot  spot'  temperature 
to  be  estimated  in  some  specified  assemblies  so  that  operational  safety 
can  be  confirmed. 

Baker  [22]  refines  methods  for  calculation  of  fuel  temperature  in 
PuOj-UC^  fuel  in  Fast  Breeder  Reactor  fuel  pins.  His  primary  concern  is 
the  calculation  of  temperature  changes  across  the  fuel-to-cladding  gap 
in  pins  with  fuel  burnups  that  range  from  60  -  10,900  MWd/MTM  (0.006  - 
1.12  at.%).  He  proposes  a  complete  set  of  heat  transfer  formulations 
that  include  sodium- to -cladding  temperature  drop,  cladding  temperature 
drop,  heat  transfer  across  the  fuel-to-cladding  gap,  and  heat  transfer 


9 


across  the  unrestructured  and  restructured  fuel  regions.  He  also  uses  a 
computer  code  (SIEX-M1)  to  calibrate  the  fuel -to- cladding  gap 
conductance  model.  Calculations  of  the  heat  transfer  across  the  fuel- 
to-cladding  gap  are  based  primarily  on  expressions  suggested  by  Ross  and 
Stoute  [23] .  The  modified  model  used  is  briefly  described  in  the 
documentation  of  the  SIEX  computer  code  [24]. 

Meek  et.  al.  [25]  apply  Kalman  filter  methodology  to  an  in-pile 
liquid-metal  fast  breeder  reactor  simulation  experiment  to  obtain 
estimates  of  the  fuel  clad  thermal  gap  conductance.  By  using  a  Kalman 
filter  (a  linear  minimum-error  variance  sequential  state  estimation 
algorithm)  an  optimal  estimate  of  the  state  vector  chosen  to 
characterize  the  experiment  is  obtained.  From  this  estimate,  the  fuel- 
clad  thermal  gap  conductance  is  calculated  as  a  function  of  time  and 
axial  position  along  the  length  of  the  fuel  pin. 

Horn  and  Fanisko  [10]  give  a  detailed  description  of  the  GAPCON 
code  and  the  different  subprograms  that  are  part  of  it.  Of  particular 
interest  to  the  present  work  is  the  subroutine  used  to  calculate 
cladding- to -coolant  heat  transfer  coefficients.  If  the  coolant  is 
water,  the  calculation  is  based  on  the  Dittus-Boelter  equation  [26],  If 
the  coolant  is  liquid  sodium,  then  the  same  expression  as  noted  by  Baker 
[22]  is  used.  The  program  is  designed  to  handle  thermal  conductivity  as 
a  function  of  temperature,  a  feature  also  available  on  the  present  work. 

Sudoh  [11]  presents  HETFEM,  a  two-dimensional  thermal  conductivity 
analysis  code  that  utilizes  a  finite  element  method  for  analyzing 
thermal  responses  of  fuel  rods  during  reflood.  He  develops  the  partial 
differential  equations  and  boundary  conditions  in  vector  notation  for 
transient  heat  conduction  in  solids  and  uses  the  method  of  weighted 


residuals  to  approximate  the  solution  to  the  problem. 

Roberts  et.  al.  [27]  present  the  Performance  Analysis  and  Design 
(PAD)  computer  code  as  an  example  of  a  computer  code  that  uses  a 
modified  version  of  the  Mikic-Todreas  model  [28]  to  calculate  gap 
conductance  after  pellet-clad  contact  has  been  established.  The  model 
operates  under  the  assumption  that,  once  contact  has  been  made,  the  gap 
conductance  is  no  longer  a  function  of  internal  pressure. 

Marr  [29]  shows  how  COBRA- 3M  was  developed  to  be  suitable  for 
analysis  of  thermal-hydraulics  in  small  pin  bundles  commonly  used  in  in¬ 
reactor  and  out-of-reactor  experiments.  It  uses  detailed  thermal 
circuits  models  for  the  fuel  pins  and  duct  walls.  The  model  handles 
temperature  dependent  properties,  nonuniform  power  distributions  and 
fuel  restructuring  and  melting.  COBRA- 3M  also  uses  finite  difference 
techniques  and  one -dimensional  (radial)  heat  flow. 

Finally,  Nguyen-Minh  and  Neuer  [30]  investigate  the  possibility  of 
using  the  modulated  heating  beam  method- -a  well-known  technique  for 
measuring  thermal  diffusivity- -for  determining  thermal  gap  conductance. 
They  find  that  if  only  radiant  heat  losses  are  present,  it  is  possible 
to  determine  thermal  gap  resistance  merely  by  measuring  the  phase  lag. 

The  preceding  summary  is  by  no  means  comprehensive.  It  is 
presented  here  to  give  the  reader  an  idea  of  the  work  being  done  in  the 
nuclear  field  pertaining  to  the  determination  of  heat  transfer 
coefficients  and  cooling  requirements  in  general.  Note  should  be  made 
that  the  emphasis  of  the  foregoing  studies  usually  centers  on  the  heat 
transfer  characteristics  of  the  gap  between  the  fuel  and  cladding.  Any 
work  done  with  melting  concerns  only  the  melting  point  of  the  fuel  since 
it  operates  at  temperatures  much  greater  than  the  cooled  cladding.  The 
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present  work,  on  the  other  hand,  is  more  concerned  with  the  cladding 
since  it  more  closely  approximates  the  "vertical  flat  plate"  of  the 
basic  problem.  But  the  data  presented  in  the  literature  is  still 
valuable  to  the  present  work  because  with  the  gap  conductance  values 
listed  and  temperature  profiles  available- -worst-case  heat  intensities 
can  be  estimated  thus  enabling  the  examination  of  conditions  leading  to 
possible  melting  of  the  cladding. 


.  ->  Mvttaciitranvt  r>  i-u  Laaci -fitittiiiiiiiit 

In  laser-machining,  Golodenko  and  Kuz'michev  [31]  examine  the 
thermal  processes  that  occur  in  metals  when  irradiated  by  powerful  laser 
pulses.  They  solve  the  heat- conduction  problem  giving  due  allowance  for 
the  motion  of  the  evaporating  boundary.  An  infinite  plane  metal  surface 
is  considered  that  has  been  subjected  to  a  uniform  surface  density  of 
radiant  energy  from  a  laser.  The  results  of  the  calculations  for  copper 
are  presented  in  the  form  of  a  series  of  graphs.  The  results  confirm 
their  expectations  that,  for  the  same  energy  on  the  pulse,  the  thickness 
of  the  evaporated  layer  is  greater  for  pulses  of  shorter  duration. 

Locke  et.  al.  [32]  examine  deep  penetration  welding  with  high-power 
C(>2  lasers  at  substantially  higher  power  levels  than  on  all  previously 
reported  experiments.  Experimental  results  are  presented  in  graph  and 
table  form.  The  measured  penetration  for  their  tests  is  in  reasonably 
good  agreement  with  established  correlations  for  vaccuum  electron  beam 
welding  [33]. 

Torvik  [16,  34]  presents  a  method  for  determining  the 


two-dimensional  transient  temperature  distribution  and  progression  of 
melting  on  disks  subjected  to  an  applied  flux  over  one  face.  He  uses  a 
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finite  element  method  that  performs  heat  balances  in  the  elements  over 
finite  time  measurements.  His  method  is  in  good  agreement  with  known 
solutions  to  two-dimensional  heat  conduction  problems  as  well  as  the 
one -dimensional  melting  problem.  As  a  result,  the  present  work  draws 
heavily  from  Torvik,  using  his  work  to  help  confirm  the  validity  of  the 
present  model  as  well  as  set  up  the  analyses.  Also,  Torvik  includes  in 
his  results  a  demonstration  that  the  time  required  for  melt-through 
approaches  that  for  the  one -dimensional  case  if  the  power  per  unit 
thickness  is  sufficiently  large.  For  future  study,  it  is  helpful  to 
note  that  phase  changes  can  also  be  incorporated  in  Torvik 's  model.  In  a 
later  numerical  study  [34],  he  produces  a  method  for  estimating  from  a 
single  curve  the  melt-through  time  for  any  arbitrary  combination  of 
flux,  spot  size,  thickness,  and  material.  Also  presented  is  melt- 
through  time  for  the  case  of  complete  retention  of  the  melt  until 
vaporization  as  well  as  predictions  for  heating  due  to  a  large  single 
pulse.  Part  of  Torvik* s  analysis  is  used  for  the  purpose  of  comparison 
in  the  present  work  since  it  most  closely  approximates  the  basic  problem 
under  study. 

Paek  and  Gagliano  [35]  develop  a  model  that  uses  a  continuous, 
distributed,  moving  heat  source  to  describe  the  temperature  profile  and 
thermal  stress  propagation  for  laser-drilled  holes  in  high-purity  fired- 
alumina  ceramic  substrate  material.  In  order  to  indicate  the  magnitude 
of  those  factors  influencing  the  potential  fracture  of  the  alumina 
material,  the  temperature  profile  and  tangential  stress  distribution  of 
the  laser-formed  hole  are  calculated.  Both  finite  and  infinite  bodies 
are  considered,  and  temperature  independent  properties  and  no  heat 
losses  through  the  bounding  surface  are  assumed.  Reflection  losses  are 
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neglected.  A  very  practical  result  of  their  work  Is  the  ability  to 
predict  the  profile  of  the  laser-drilled  hole.  The  theory  they  develop 
indicates  shorter  pulse  lengths  generally  give  lower  stresses  and  that 
the  desirable  condition  for  hole  drilling  is  to  maintain  a  sufficiently 
high  beam  intensity  for  a  minimum  time  in  order  to  remove  by  sublimation 
a  given  mass  of  material. 

Prokhorov  et.  al.  [36]  present  experimental  and  theoretical  results 

on  metal  evaporation  under  powerful  optical  radiation  (normal  laser 
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pulses  with  light  intensity  of  approximately  10  W/m  and  spot  size  on 
the  target  of  about  0.01  m) .  They  also  review  some  factors  that  limit 
the  rise  of  surface  temperature  (i.e.  instability  of  superheated 
liquids,  intense  boiling  of  liquid,  the  disappearance  of  metal 
conductivity).  Furthermore,  they  point  out  that  experimental  data  is 
already  in  the  literature  directly  supporting  the  existence  of  a  sharp 
drop  in  reflectivity  with  increased  intensity  of  light  incident  on  a 
metal  surface  [37].  They  propose  a  method  for  calculating  metal  surface 
temperature  versus  incident  light  intensity  by  using  an  approximate 
solution  to  the  Clapeyron-Clausius  equation.  Results  of  investigating 
vapor  dynamics  of  metal  evaporation  under  powerful  millisecond  optical 
radiation  are  also  reported. 

Jones  and  Vang  [38]  examine  the  effect  of  temperature  as  a  function 
of  time  in  laser  welding  of  similar  and/or  dissimilar  materials  using 
experiment  and  finite  element  analyses.  They  detail  the  laser  operating 
parameters  as  well  as  the  finite  element  method  and  temperature 
dependent  properties.  Their  experimental  results  show  that,  with  proper 
operating  parameters,  good  tensile  strength  and  low  electrical 
rfsistance  can  be  obtained  in  copper-copper  welds.  With  the  addition  of 
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moderate  contact  pressure,  copper -aluminum  welds  can  result  in  good 
tensile  and  impact  strength  as  well  as  good  electrical  resistance 
stability  after  thermal  cycling.  For  the  present  work,  their  use  of 
finite  element  analysis  is  the  most  noteworthy. 

Blottner  [39]  presents  a  relatively  simple  numerical  technique  for 
predicting  the  melt  depth  for  pulsed  laser  welding.  He  approximates  the 
two-dimensional  aspects  of  the  problem  by  using  an  effective  area  in  the 
one -dimensional  energy  equation.  He  uses  the  enthalpy  method  with  an 
implicit  numerical  scheme  to  provide  a  rapid  solution  procedure. 

Herziger  [16]  examines  the  strong  coupling  between  laser  and  target 
due  to  optical  feedback  from  the  target  that  results  in  changing  the 
system  parameters  of  the  laser  tool.  He  studies  the  processing 
parameters;  factors  affecting  the  absorption  of  laser  radiation;  and  the 
effects  of  laser  induced  plasma.  He  concludes  that  the  application  of 
absportive  layers  or  an  appropriate  choice  of  laser  parameters  for  a 
poor  optical  absorptivity  material  can  increase  its  absorptivity  to 
nearly  unity.  For  the  present  work,  this  conclusion  is  used  as  a  basis 
for  neglecting  plasma  effects.  He  further  concludes  that  the  quality  of 
laser  machining  depends  strongly  on  the  proper  adjustments  and  control 
of  the  laser  radiation. 

Decker  et.  al.  [40]  present  some  primary  steps  in  developing  a 
cutting  model.  They  conclude  a  cutting  model  should  be  as  simple  as 
possible  and  should  be  able  to:  give  an  estimation  of  the  cutting 
speed;  show  how  the  quality  of  the  cutting  surface  depends  on  the 
process  parameters;  and  predict  the  metallurgical  impact  on  the 
workpiece.  (For  the  present  work,  these  attributes  of  a  good  cutting 
model  were  the  deciding  factors  in  choosing  the  work  of  Torvik  [16]  as 
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the  basis  of  comparison  since  he  not  only  uses  finite  element  analysis 
but  his  model  addresses  cutting  speed  and  the  metallurgical  impact.) 
Decker  et.  al.  also  list  the  three  different  versions  of  laser  cutting 
as  well  as  the  parameters  controlling  the  process  (i.e.  laser  beam, 
assistant  gas,  material,  slag).  The  need  to  optimize  the  parameters  of 
cutting  speed,  width,  and  cutting  surface  is  presented  in  addition  to 
details  on  each  parameter.  The  cutting  process  is  thoroughly  examined 
and  the  report  is  concluded  with  a  detailed  cutting  model  that  considers 
the  dissipation  and  production  of  heat  energy. 

As  with  nuclear  fuel,  the  above  listing  is  by  no  means  exhaustive. 
The  listing's  primary  purpose  is  to  give  the  reader  a  general  overview 
of  work  already  accomplished  in  laser-machining  that  applies  to  the 
present  work  and  to  some  of  its  possible  applications.  Additionally, 
the  listing  presents  a  brief  introduction  to  the  fundamental  nature  of 
the  basic  problem  of  localized  high  intensity  heat  sources. 


Chapter  3 
MODEL  DEVELOPMENT 


3.1  Preliminaries 

The  first  decision  that  had  to  be  made  was  to  choose  which  computer 
system  to  use  for  development  of  the  model.  The  logical  choice  was  the 
University  of  Kentucky  Engineering  Computer  Center's  Harris  800  Computer 
System.  This  system  had  just  recently  been  donated  to  the  college  by 
the  Harris  Corporation  and  as  a  result  the  computer  time  was  'free'. 

With  unlimited  computer  time  available  coupled  with  the  fact  that  this 
model  was  to  provide  only  the  foundation  for  more  in-depth  examination 
of  the  basic  problem,  the  next  decision  was  to  proceed  with  a  three- 
dimensional  model.  Most  of  the  literature  deals  with  one-  and  two- 
dimensional  models  due  to  the  normal  high  expense  of  computer  time  and 
the  need  for  expeditious  solutions  to  engineering  problems.  With 
neither  of  these  constraints  present,  a  three-dimensional  model  would 
allow  the  greatest  accuracy  as  well  as  the  greatest  versatility  since  it 
could  examine  cases  where  the  one-  and  two-dimensional  approximations 
would  not  apply. 

With  a  computer  system  and  dimension  chosen,  a  methodology  had  to 
be  determined.  The  most  extensively  used  method  is  probably  the  finite 
difference  method.  Its  use  in  one-,  two-,  and  three-dimensional  models 
is  well  documented  [41-45],  and  involves  the  use  of  discretization  to 
develop  algebraic  "difference"  equations  from  the  governing  equations. 
Advantages  of  this  numerical  method  are  its  familiarity  (most  computer 
courses  cover  finite  differences  before  other  numerical  methods)  and  its 
relative  simplicity.  A  disadvantage  of  this  numerical  method,  as  with 
any  method  designed  to  solely  calculate  numerical  values  for 
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temperatures,  is  that  temperature  gradients  cannot  be  determined 
accurately.  Thus,  thermal  stresses,  which  depend  directly  on  thermal 
gradients,  cannot  be  easily  calculated.  Also,  phase  change  and  material 
inhomogeneities  are  difficult  to  model. 

Alternatively,  the  finite  element  method  offers  some  different 
features.  It  too  has  been  well  documented  in  one-,  two-,  and  three- 
dimensional  use  [11,  16,  34,  38,  46  -  49],  and  in  addition  is 
advantageous  in  dealing  with  inhomogeneous  materials  and  phase  changes. 
Additional  advantages  of  the  finite  element  method  over  other  numerical 
approaches  are  listed  by  Wilson  and  Nickell  [46].  For  the  finite 
element  method,  the  body  under  study  is  replaced  by  a  system  of  elements 
with  each  element  being  defined  by  a  pattern  of  nodes.  For  heat 
transfer  studies,  an  approximate  solution  for  each  element  in  the 
temperature  field  is  assumed  and  heat  flux  equilibrium  equations  are 
developed  within  the  finite  element  system  at  a  discrete  number  of 
nodes.  The  algebraic  relationships  between  the  nodal  temperatures  are 
then  solved.  For  transient  problems,  the  system  of  equations  must  be 
solved  for  each  time  step  which  requires  the  inversion  of  a  large  matrix 
for  each  time  Increment.  This  can  lead  to  long  computer  running  times 
but  the  flexibility  of  the  method  more  than  makes  up  for  this 
disadvantage . 

Weighing  the  above  considerations  made  the  choice  of  the  finite 
element  method  the  likely  result,  but  the  decisive  factor  was  the 
availability  of  the  finite  element  program  ANSYS  by  Swanson  Analysis 
Systems,  Inc.  on  the  Harris  computer.  Thus,  the  combination  of  a  ready¬ 
made  finite  element  system  and  unlimited  computer  time  made  finalizing 
the  preliminary  decisions  a  rapid  procedure. 


3.2  ANSYS 

The  ANSYS  computer  program  is  a  large-scale,  general  purpose 
program  for  the  solution  of  several  varieties  of  engineering  analyses. 
The  program  employs  the  matrix  displacement  method  of  analysis  based  on 
the  finite  element  idealization.  It  uses  a  wave-front  (or  frontal) 
solution  procedure  [50,  51]  for  the  system  of  simultaneous  linear 
equations  developed  from  the  set  of  finite  elements.  The  wave  front  at 
a  point  is  defined  as  the  number  of  active  equations  after  any  element 
has  been  processed  during  the  solution  procedure.  The  program  is 
designed  to  handle  up  to  2000  degrees  of  freedom  normally,  however  the 
version  used  by  universities  is  a  special  educational  version  of  ANSYS 
and  is  limited  to  a  200  degree  of  freedom  wave  front.  This  and  the 
limited  memory  size  for  node  coordinate  storage  are  the  only  differences 
between  the  two  versions,  but  they  allow  Swanson  Analysis  Systems,  Inc. 
to  lease  ANSYS  to  the  schools  at  a  greatly  reduced  cost. 

Though  ANSYS  is  mainly  a  structural  analysis  program,  it  handles 
heat  transfer  problems  equally  well.  Its  thermal  analysis  module  is 
used  to  solve  for  steady-state  or  transient  temperature  distributions  in 
a  body.  Conduction,  convection,  radiation,  and  internal  heat  generation 
can  all  be  incorporated  in  the  model  along  with  heat  transport  effects 
due  to  flowing  fluids.  Material  properties  may  be  either  constant  or 
temperature  dependent,  and  output  temperatures  may  be  stored  and  used 
for  structural  analyses  (a  feature  that  will  be  utilized  in  future 
work).  Automatic  convergence  criteria  are  available  for  steady-state, 
nonlinear  solutions  while  a  time -step  optimization  procedure  is 
available  for  transient  solutions.  For  the  present  work,  both  of  these 


features  are  demonstrated.  The  latter  one  is  used  to  verify  the  model, 
while  the  former  is  used  for  the  major  part  of  the  analysis. 

The  basic  equation  for  the  thermal  analysis  in  ANSYS  is  Poisson's 
equation  with  temperature  being  the  primary  unknown.  The  basic  equation 
is  written  as: 

[C]{f]  +  [K]{T)  -  {Q} 

where : 

[C]  -  the  specific  heat  matrix 

[K]  -  the  thermal  conductivity  matrix  (including  equivalent  face 
convection  terms) 

(T)  -  the  nodal  temperature  vector 

{ Q )  -  the  heat  flow  rate  vector  (including  applied  heat  flow, 
internal  heat  generation,  and  convection) 

(Note:  For  the  steady-state  case,  which  is  the  main  concern  of  the 
present  work,  the  (f)  term  is  zero  and  the  equation  reduces  to:  [K] (T ) 

-  { Q ) . )  The  basic  equation  is  solved  by  an  implicit  direct  integration 
scheme  based  on  a  modified  Houbolt  method.  This  method  uses  a  quadratic 
temperature  function  (where  T  is  the  temperature  at  time  t)  over  the 
region  being  examined  of  the  form: 

2 

T  -  a  +  bt  +  ct 

The  temperature  function  is  differentiated  and  substituted  into  the 
basic  equation  to  yield  an  equation  with  three  unknowns:  a,b,  and  c. 

If  At  is  the  time  step  between  iterations,  a  set  of  three  equations  may 
be  defined  at  t,  t  -  At,  and  t  -  2At.  These  equations  are  then  solved 
simultaneously  to  give  the  integration  equation: 
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M1 

At 


[C]  [K] 


<V  -  (Qt>  +  f([C].(Tt.At},(Tt.2At 


)) 


(where:  is  a  constant  and  Qt  is  the  heat  flow  rate  vector  at  time  t.) 

This  equation  may  be  solved  for  T  since  solutions  at  previous  times  are 
known.  For  more  detailed  information  on  the  working  of  ANSYS,  the 
reader  is  referred  to  the  appropriate  manuals:  DeSalvo  [52]  and  Kohnke 
[53]. 


3.3  Verification  of  ANSYS 

Once  ANSYS  was  chosen  as  the  finite  element  program  to  be  used  in 
the  model,  it  had  to  be  tested  and  verified.  The  work  selected  for 
comparison  was  that  of  Torvik  [16]  since  his  finite  element  model  was 
designed  t-o  handle  phase  change  along  with  temperature  distribution. 
Also,  he  gives  a  detailed  comparison  with  known  analytical  solutions. 

The  third  example  in  his  work  is  the  case  of  the  two-dimensional 
temperature  distribution  resulting  from  a  steady  axi- symmetric  flux  of  a 
Gaussian  distribution  acting  on  a  thin  sheet.  The  Gaussian  distribution 
is  of  the  form: 

F(r)  -  F  x  exp(-r2/2o2) 

Torvik  compares  his  solutions  with  the  temperature  profiles  predicted  by 
a  computer  program  written  at  the  Air  Force  Weapons  Laboratory  which 
evaluates  the  analytical  Fourier  series  solution  given  by  Olcer  [54]. 

A  specific  case  Torvik  considers  is  that  of  a  titanium  sheet, 

0.0004  m  thick  and  0.05  m  in  diameter,  initially  at  300*K,  and  having 
the  following  thermal  properties: 
p  -  4430  kg/m^ 


k  -  14.5  W/(m*C) 
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cp  -  770  Joules/(kg#C) 

T  -  1900°K 
m 

He  assumes  the  disk  is  subjected  to  a  beam  of  peak  intensity  of  2.0 
7  2 

x  10  V/m  with  a  -  0.0025  m.  Torvik  uses  his  model  to  examine  the 
results  of  dividing  the  slab  into:  a)  10  layers  and  40  annular  rings 
using  a  time  step  of  80  psec,  and  b)  20  layers  and  20  annular  rings 
using  a  time  step  of  20  fisec.  Both  cases  yielded  temperature  profiles 
that  differed  by  less  than  1%  at  a  depth  of  2.0  x  10  ^  m  for  all  r  and 
at  times  of  t  -  0.03  and  0.10  sec.  Also,  both  cases  were  in  excellent 
agreement  with  the  Fourier  series  solution  program,  especially  within 
the  beam  radius  (2o  -  0.005  m) .  The  exact  solution  indicated  that 
melting  begins  at  the  front  surface  at  t  -  0.1037  sec.  Torvik' s  finite 
element  solutions  gave  0.1042  and  0.1064  sec.,  respectively. 

With  these  data  for  comparison,  the  ANSYS  model  was  set  up  to 
establish  its  accuracy.  The  Torvik  problem  was  changed  only  slightly  to 
facilitate  later  analysis  of  results  using  flat  plate  correlations. 

Thus,  a  0.05  m  square  plate  was  substituted  for  the  0.05  m  radius  disk. 
All  other  parameters  remained  the  same.  Since  the  present  work  is 
concerned  only  with  preventing  the  plate  from  starting  to  melt,  the 
crucial  comparison  became  the  exact  solution's  time  of  t  -  0.1037  sec. 
(when  melting  begins  at  the  front  surface) . 

Because  the  educational  version  of  ANSYS  is  limited  to  only  a  200 
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Torvik's  case  b) ,  the  equivalent  of  20  annular  rings  was  set  up- -namely, 

a  20  by  20  element  arrangement  over  the  quarter  of  the  plate.  This 

-3 

resulted  in  400  square  elements  that  were  1.25  x  10  m  per  side 
covering  the  quarter  of  the  plate.  However,  since  the  educational 
version  of  ANSYS  could  not  accommodate  10  layers  (much  less  20  layers) 
due  to  the  wave-front  limitation,  4  layers  had  to  be  selected.  The 
Gaussian  distribution  of  the  Torvik  problem  was  handled  by  taking  the 
radial  distance  of  the  center  of  the  face  of  each  element  within  the 
beam  radius,  substituting  this  distance  as  r  into  the  Gaussian 
distribution  equation,  and  using  the  resultant  intensity  value  as  the 
average  heat  transfer  coefficient  over  that  element  face.  Since  square 
elements  were  being  used  to  approximate  a  circular  beam  area,  an  area 
3.45%  larger  had  to  be  used.  The  resultant  effect  was  expected  to  lower 
the  melting  initiation  time  only  a  slight  amount  (at  most  a  proportional 
percentage)  due  to  the  fact  this  extra  area  was  in  the  edges  of  the 
affected  area  on  the  plate- -the  region  of  lowest  intensity  because  of 
the  beam's  Gaussian  distribution.  For  comparison,  the  model  was  run 
under  two  additional,  different  conditions:  a)  the  central-most  element 
within  the  beam  radius  (in  this  case  the  lower-left  corner  element)  was 

given  a  convective  boundary  condition  value  equal  to  the  peak  value  of  2 
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x  10  Vatts/m  (this  amounted  to  an  increase  of  1.2  x  10  Watts/m  ) 

which  was  estimated  to  give  a  better  approximation  due  to  the  nature  of 

the  Gaussian  distribution;  and  b)  the  entire  beam  area  was  given  a 

convective  boundary  condition  equal  to  the  peak  value  in  order  to  note 

the  error  incurred  by  a  uniform  peak  loading. 

The  results  are  presented  in  Table  1  and  indicate  the  ANSYS  model 
is  quite  acceptable,  especially  when  using  the  better  approximation  to 
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TABLE  1 

Test  Runs  to  Establish  Accuracy  of  ANSYS  Model1 


Tine  at 
Initiation 

of  Melting  %  error2 

(sec) 


Temperature 
at  t  -  0.1037 
[exact  solution]  %  error 
(eK) 


1)  1.88xl07  W/m2 
on  central 
elementc4 

0.1095 

+  5.60% 

1829° 

-  3.74% 

2)  2 . OxlO7  W/m2 
on  central 
element4 

0.1022 

-  1.45% 

1919° 

+  1.00% 

3)  2. OxlO7  W/m2 
uniform  loading 
on  all  elements 
in  affected 

area 

0.0970 

-  6.46% 

1998° 

+  5.16% 

1  Note:  Cases  1)  and  2)  above  use  the  values  noted  for  the  central-most 

element,  and  the  average  values  based  on  radial  distance  to  the 
center  of  the  face  of  the  element  for  the  other  elements  within 
the  beam  radius . 

2  Note:  Percent  error  is  based  on  comparison  to  exact  solution  of  t  - 

0.1037  s  for  initiation  of  melting  time. 

3  Note:  Percent  error  is  based  on  comparison  to  the  melting  point 

temperature  of  titanium  used  by  Torvik  [16]  :  1900°K. 

4  Note:  Since  the  symmetry  of  the  problem  is  being  invoked,  only  the 

upper-right  quarter  of  the  plate  is  being  examined.  Thus,  the 
central  element  is  in  actuality  the  bottom- left  corner  element 
on  the  portion  of  the  plate  being  examined.  Also,  this  central 
element  represents  1/4  of  the  central-most  actual  area  of  the 
incident  Gaussian  laser  Learn. 
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the  Gaussian  distribution  (peak  intensity  for  central-most  element). 

Since  it  was  expected  that  the  3.45%  area  increase  would  yield  a 

slightly  lower  initiation  of  melting  time,  the  results  for  the  case  of  2 
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x  10  Watts/m  are  especially  encouraging  and  indicate  that  the  ANSYS 
model  is  within  about  1%  of  the  analytical  solution.  Also  noteworthy  is 
the  fact  that  practically  no  accuracy  was  lost  due  to  using  4  layers  of 
elements  as  opposed  to  Torvik’s  10  and  20  layer  arrangements.  Of  course 
this  result  can  only  be  applied  with  certainty  to  the  case  of  0.0004  m 
thick  plates  without  more  testing.  However,  since  the  present  work  is 
limited  to  this  thickness  (except  for  an  example  case  using  aluminum), 
this  finding  only  further  confirms  the  applicability  and  accuracy  of 
ANSYS  to  the  portion  of  the  basic  problem  under  consideration.  For 
verification  purposes,  a  copy  of  one  of  the  computer  programs  used  to 
test  the  accuracy  of  ANSYS  is  submitted  in  Appendix  A. 

3.4  Defining  Model  Parameters 

With  ANSYS  confirmed  as  appropriate  for  examining  the  basic 
problem,  the  task  undertaken  next  was  defining  the  limits  of  the 
problem.  Since  ANSYS  had  just  been  tested  using  Torvik's  third  example, 
his  parameters  were  retained.  Thus,  the  plate  under  consideration  took 
on  the  dimensions  of  0.05  m  by  0.05  m  with  a  depth  of  0.0004  m.  For  the 
materials  to  be  studied,  titanium  was  an  obvious  choice  since  it  had 
just  been  tested  and  since  it  has  many  practical  applications  in  the 
nuclear  and  laser-machining  fields.  For  comparison,  two  other  widely 
used  materials,  aluminum  and  stainless  steel  (AISI  304),  were  also 
chosen. 

Aluminum  was  selected  in  order  to  examine  a  metal  with  considerably 
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different  properties  than  titanium.  Aluminum's  melting  point  is  nearly 
1000#K  below  that  of  titanium  while  its  conductivity  is  approximately  an 
order  of  magnitude  greater.  Stainless  steel,  on  the  other  hand,  was 
selected  because  its  properties  fell  in  between  those  of  titanium  and 
aluminum- -although  they  were  more  on  the  order  of  titanium. 

With  the  dimensions  and  materials  of  the  plate  established,  the 

7  2 

heating  conditions  were  next.  Since  2  x  10  W/m  had  just  been  tested, 

it  was  chosen  as  the  base  value  of  intensity.  Since  Table  1  indicated  a 

+5.16%  error  in  temperature  was  obtained  when  using  a  uniform  loading  of 

peak  intensity,  uniform  loading  was  chosen  to  insure  a  safety  factor  as 

well  as  determine  an  upper  bound.  Torvik's  beam  radius  (0.005  m)  was 

retained  as  a  base  value  so  that  the  initial  affected  area  for  the 

quarter  of  the  plate  consisted  of  13  elements  [total  area  -  4  x  13  x 
-32  -52 

(1.25  x  10  m)  -  8.125  x  10  m  ] .  In  order  to  examine  the  effects  of 
area  and  intensity  on  cooling  requirements,  the  intensity  was  chosen  to 
vary  from  1  x  10^  W/m2  to  6  x  102  W/m2  in  steps  of  1  x  10^  W/m2. 
Affected  area  was  also  varied  from  approximately  30.8%  of  the  original 
area  (2.5  x  10*5  m2)  to  400%  (3.25  x  10*4  m2)--500%  (4.0625  x  10'4  m2) 
for  aluminum- -with  the  combined  results  being  presented  in  tabular  and 
graphical  form. 

For  the  cooling  part  of  the  problem,  liquid  sodium,  air,  and  water 
were  chosen  as  the  coolants.  Liquid  sodium  was  selected  due  to  its  low 
Prandtl  number  and  air  was  selected  as  a  representative  gas.  It  was 
desired  to  also  choose  a  fluid  of  high  Prandtl  number,  but  at  the 
average  temperatures  the  model  deals  with  (600-950°K) ,  all  high  Prandtl 
number  fluids  the  author  could  identify  decomposed.  Thus,  liquid  water 
at  high  temperatures  and  pressures  was  chosen  in  order  to  examine  at 
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least  one  case  of  relatively  high  Prandtl  number  (Pr  -  12.0). 

To  further  establish  upper-bound  conditions,  all  cooling  effects  on 
the  surface  where  the  heat  source  was  located  were  neglected  as  well  as 
any  losses  out  the  plate  edges.  Also,  the  model  examined  only  the 
steady-state  condition  so  as  to  provide  an  upper-bound  temperature 
distribution.  (Note:  For  those  cases  where  the  steady-state  is  not  an 
upper  bound,  the  intensity  of  the  heat  source  will  be  so  great  that 
melting  will  occur  regardless  of  cooling  conditions.  Torvik's  first 
example  [16]  demonstrates  this  by  showing  that  a  finite  slab  can  be 
considered  as  being  of  infinite  thickness  for  very  short  initial  time 
periods- -0.06  seconds  for  his  example.)  The  rationale  behind  the  above 
assumptions  was  that  if  cooling  conditions  could  be  established  that 
would  prevent  melting  under  these  upper-bound  conditions,  then  these 
cooling  conditions  should  offer  a  margin  of  safety  under  actual 
conditions  as  well  as  reliable  minimum  cooling  requirements. 

Thus,  with  all  these  factors  incorporated,  the  final  product  of  the 
model  becomes  a  heat  transfer  coefficient  that  will  prevent  melting 
under  upper-bound  conditions.  This  heat  transfer  coefficient  is  then 
used  along  with  local  Nusselt  number  correlations  for  a  vertical  flat 
plate  to  determine  if  cooling  is  possible  under  free  or  forced 
convection  and,  if  cooling  is  indeed  possible,  under  what  specific 
conditions.  A  copy  of  a  sample  program  used  in  determining  these 
conditions  is  submitted  in  Appendix  B.  The  results  of  the  model's 
computations  are  presented  and  discussed  in  the  next  chapter. 


Chapter  4 

NUMERICAL  SOLUTIONS  AND  CORRELATION  RESULTS 


4.1  Presentation  of  Numerical  Solutions 

In  order  to  obtain  the  desired  solutions,  the  author  had  to 
function  as  an  optimizer  for  the  ANSYS  program- -supplying  estimates  of 
the  appropriate  heat  transfer  coefficient  to  the  program  and  then 
altering  the  value  after  the  program  had  returned  an  appropriate 
temperature  distribution.  (Swanson  Analysis  Systems,  Inc.,  just  recently 
added  an  optimizer  module  to  ANSYS.  Unfortunately,  It  was  not  available 
during  the  execution  of  this  work. )  Since  the  average  run  took  about  40 
minutes  of  CPU  time  (and  up  to  three  hours  of  connect  time),  the  process 
proved  to  be  extremely  time  consuming.  Fortunately,  the  author  was 
given  four  different  computer  accounts  and  this  expedited  matters 
somewhat.  Even  so,  the  above  time  constraints  forced  some  interpolation 
of  values  that  the  author  would  have  preferred  to  forgo,  however  these 
are  noted  where  they  occur  and  are  estimated  to  be  well  under  1%  in 
error.  Summaries  of  the  results  are  presented  below. 

4-1-1  Heat  Transfer  Coefficient  as  a  Function  of  Heat  Source  Intensity 

Figure  1  is  a  combined  graph  of  the  response  of  all  three  metals  to 
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the  range  of  intensities  (1  -  6  x  10  W/m  )  under  study.  The  initial 
response  is  as  intuitively  expected,  with  the  metals  having  the  highest 
conductivities  requiring  the  most  cooling.  The  approximate  linearity  of 
the  aluminum  response  can  be  partly  explained  by  its  high  conductivity 
(at  least  an  order  of  magnitude  greater  than  titanium  and  stainless 
steel)  and  partly  explained  by  the  fact  its  conductivity  is  decreasing 
with  temperature.  The  similarity  of  the  stainless  steel  and  titanium 


ALUMINUM  (MP  -  933  K) 
STAINLESS  STEEL  (MP  -  1670  K) 


Minimum  heat  transfer  coefficient  [h(x)J  plotted  as  a  function  of 
localized  heat  source  intensity  for  aluminum,  stainless  steel,  an< 
titanium  plates,  0.0004  m  thick. 


responses  are  expected  since  their  conductivities  are  approximately  the 

same  and  also  increase  with  temperature. 

To  fully  understand  Figure  1,  especially  the  aluminum  response 

curve's  crossing  of  the  titanium  and  stainless  steel  response  curves 
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between  5  -  6  x  10  W/m  ,  an  examination  of  the  basic  governing 
equations  is  necessary.  The  heat  flux  (q^)  through  the  plate  is 
approximately  defined  by  a  simplified  form  of  Fourier's  law  to  be: 

Ox  -  KTj.  -  VA 

where:  k  is  the  conductivity;  T_  is  the  temperature  on  the  front  face 

r 

of  the  plate  (the  side  with  the  heat  source);  Tfi  is  the  temperature  on 
the  back  face  of  the  plate  (the  side  with  the  coolant);  and  L  is  the 
thickness  of  the  plate. 

For  the  convection  process  on  the  back  of  the  plate,  Newton's  law 
of  cooling  defines  the  heat  flux  from  the  plate  (in  this  case,  q^)  to 
be : 

Ox  -  MTb  -  V 

where:  h  is  the  heat  transfer  coefficient,  Tfi  is  the  temperature  of  the 
back  face  of  the  plate;  and  T^  is  the  fluid  free  stream  temperature. 

Combining  and  rearranging  the  above  two  relations  to  obtain  an 
expression  for  h  yields: 

h  -  k(Tp  -  Tb)/L(Tb-TJ 

For  the  case  of  each  metal  in  the  present  work:  T^,,  L  and  T  are  all 
constant.  The  conductivity  (k)  varies  with  temperature  through  the 
plate  for  all  three  metals,  but  since  k  is  so  large  for  aluminum,  the 
temperature  difference  through  the  aluminum  plate  is  relatively  constant 
(varying  from  33*  to  about  100°)  for  all  intensities.  This  explains  the 
resulting  near -linear  behavior. 


On  the  other  hand,  for  the  case  of  titanium  and  stainless  steel, 

not  only  does  k  increase  with  temperature  (and  thus  intensity  of  heat 

source),  but  TD  decreases  dramatically  with  increase  of  intensity.  For 

example,  the  temperatures  on  the  back  face  of  titanium  and  stainless 

steel  are  approximately  1600*K  and  1400°K,  respectively,  for  the  2  x  10^ 
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W/m  intensity.  For  6  x  10  W/m  ,  the  temperatures  have  dropped  to 

762. 5*K  and  665*K  respectively.  Thus,  not  only  is  the  denominator 

decreased,  but  the  numerator  is  markedly  increased  too,  so  that  the 

value  of  h  increases  almost  exponentially.  This  behavior  can  be 

attributed  to  the  relatively  poor  conductivity  of  stainless  steel  and 

titanium  as  compared  to  aluminum.  Thus,  for  stainless  steel  and 

titanium  to  keep  from  melting  as  the  intensity  increases,  a  large 

temperature  gradient  must  be  created  by  almost  exponential  increases  of 

heat  transfer  by  convection.  Therefore,  once  examined  in  detail,  the 

crossing  of  the  response  curves  is  expected. 

7  2 

Note  must  also  be  made  of  why  6  x  10  W/m  was  determined  as  the 

maximum  intensity  to  be  examined.  Since  the  present  work  is  concerned 

7  2 

with  only  free  and  forced  convection  cooling,  exceeding  6  x  10  W/m 
would  definitely  produce  heat  transfer  coefficients  outside  of  those 
cooling  regimes.  According  to  Incropera  and  DeWitt  (55),  typical  values 
of  heat  transfer  coefficient  for  forced  convection  using  liquids  range 
from  50  -  20,000  W/m^"K  while  boiling  values  range  up  to  100,000  W/m^°K. 

Since  the  model  produces  heat  transfer  coefficients  in  excess  of  100,000 

2  7  2 

W/m  °K  for  a  heat  intensity  of  6  X  10  W/m  this  intensity  is  easily  an 

upper  bound  for  the  present  work. 

A  listing  of  the  model  generated  data  for  Figure  1  appears  in 

Appendix  C.  Conductivity  values  are  confirmed  in  Incropera  and  Dewitt 


[55]  and  Touloukian  and  Ho  [56]. 


4.1.2  Heat  Transfer  Coefficient  as  a  Function  of  Affected  Area  at 
Various  Intensities 

Figures  2-4  show  the  combined  relations  between  heat  transfer 

coefficient  as  a  function  of  affected  area  at  various  intensities 

7  2  7  2 

ranging  from  1  -  6  x  10  W/m  .  The  2  x  10  W/m  intensity  line  is  the 

base  line  in  each  graph  since  it  has  the  most  complete  set  of  data  for 

7  2 

each  material.  Data  points  for  all  2  x  10  W/m  lines  are  plotted  as 
triangles.  Squares  mark  the  value  of  other  data  points  that  were 
actually  computed  by  the  model  either  by  exact  calculation  or 
interpolation/extrapolation  of  model  generated  values.  The  data  from 
Figure  1  is  also  included  and  plotted  as  squares. 

All  graphs  behave  similar- -with  a  fairly  rapid  rise  as  the  affected 
area  first  begins  to  grow  to  an  almost  just  as  rapid  leveling-off  (after 
a  certain  "critical"  area)  beyond  which  increase  of  area  has  relatively 
little  effect  on  heat  transfer  coefficient.  The  apparent  "critical" 
area  is  reached  earlier  with  lower  conductivity  values,  probably  because 
the  conditions  that  approximate  the  case  of  a  constant  heat  flux  over 
the  entire  surface  are  arrived  at  sooner.  This  is  probably  the  reason 
for  the  fairly  rapid  leveling-off  behavior  too  since,  as  the  area 
affected  by  the  heat  source  increases ,  the  problem  begins  to  locally 
take  on  the  appearance  of  a  wall  over  which  its  entire  surface  is 
subjected  to  a  constant  heat  flux.  From  a  local  perspective  with 
respect  to  the  plate,  the  case  of  constant  heat  flux  is  exactly  what  is 
transpiring,  and  so  a  nearly  constant  value  of  heat  transfer  coefficient 
is  to  be  expected. 

Another  interesting  trend  shown  on  the  charts  is  found  in  the 
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graphing  of  the  study's  intensity  range  (1  -  6  x  10  W/m  )  values.  In 
the  case  of  stainless  steel  and  titanium,  with  their  relatively  low  heat 
conductivities  that  increase  with  temperature,  an  ever-increasing 
spacing  is  noted  as  intensity  increases.  For  aluminum,  with  its 
relatively  high  conductivity  that  decreases  with  temperature,  the 
spacing  is  approximately  equal.  In  light  of  Figure  1,  this  graphical 
result  was  expected. 

In  addition  to  these  values,  the  model  calculated  some  key  points 
that  are  also  plotted  as  squares  on  the  graph  in  order  to  predict  the 
shapes  of  the  rest  of  the  curves .  Some  of  the  values  are  generated  by 
the  model  exactly  while  others  are  interpolated/extrapolated  from  model 
data.  Data  points  used  are  tabulated  in  Appendix  D  with  interpolations 
and  extrapolations  noted. 

4.1.3.  Heat  Transfer.  Cpefficlent  .gs.B.  Functipn  pf  Ihlgknep?  in 
AlmiiinwR 

During  the  course  of  the  present  work,  it  was  discovered  that 
because  of  aluminum's  high  conductivity,  the  4-layer  element  arrangement 
could  be  replaced  by  just  one  layer  with  a  change  of  only  about  0.1'K. 

In  fact,  it  was  learned  that  the  single  element's  depth  could  be 
extended  to  0.0008  m  with  less  than  a  0.2*K  change.  This  discovery  not 
only  cut  down  the  run  time  of  the  aluminum  jobs  to  about  25%  of  the 
original  time,  but  also  allowed  the  examination  of  aluminum  plates  up  to 
0.0032  m  thick  while  maintaining  essentially  the  same  accuracy  (within 
0.02%).  The  graph  showing  the  findings  is  presented  in  Figure  5,  while 
tabulated  data  is  submitted  in  Appendix  E.  The  best-fit  curve  presented 
is  a  regression  polynomial  of  order  2.  The  one  data  point  denoted  by  a 
square  is  the  theoretical  value  of  the  heat  transfer  coefficient  for  the 


case  of  a  plate  with  ’zero'  thickness  (i.e.  100%  heat  conduction).  The 

extrapolated  value  of  the  best-fit  curve  with  the  y-axis  yields  a  value 
4  2 

of  about  3.38  x  10  W/m  °K--or  less  than  7.0%  high.  The  difference  can 
be  explained  by  the  fact  that  the  idealization  of  using  the  average  of 
the  plate  temperature  and  the  fluid  free  stream  temperature  to  determine 
coolant  properties  will  generally  yield  local  temperatures  lower  than  in 
the  actual  case.  Still,  the  curve  behaves  as  expected,  with  the  heat 
transfer  coefficient  decreasing  with  increased  plate  thickness. 

4.1.4  Heat  Transfer  Coefficient  Versus  Conductivity 

A  bar  graph  of  the  minimum  heat  transfer  coefficients  for  the 

various  metals  versus  each  metal's  conductivity  is  presented  in  Figure 

6.  The  values  are  those  which  result  from  the  base  intensity  of  2  x  107 
2 

W/m  .  Only  three  model -generated  data  points  are  available  at  this 
intensity  since  only  three  materials  were  studied.  Further  study  of 
more  materials  is  required  to  establish  possible  trends.  The  data  for 
this  graph  is  submitted  in  Appendix  F. 

4.2  Nusselt  Number  Correlations  for  Free  and  Forced  Convection 

Given  the  heat  transfer  coefficients  needed  to  prevent  melting,  it 
was  now  necessary  to  translate  them  into  suggested  cooling  conditions. 

To  do  this,  Nusselt  number  correlations  were  required  since  the  Nusselt 
number  provides  a  measure  of  the  convection  heat  transfer  occurring  at 
the  surface  of  the  plate.  Thus,  Nusselt  number  correlations  were  needed 
for  a  vertical  flat  plate  under  both  free  and  forced  convection 
conditions.  Even  though  according  to  Incropera  and  DeWitt  [55]  the  heat 
transfer  coefficients  produced  by  the  model  were  in  the  upper  range  of 
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Figure  6.  Minimum  heat  transfer  coefficient  versus  each  metal's  conductivity 
(titanium,  stainless  steel,  and  aluminum)  for  the  case  of  heat 
intensity  -  2  x  107  W/m2 . 
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typical  forced  convection  values  for  liquids,  the  free  convection 
correlations  were  still  investigated  mainly  for  the  purpose  of  examining 
the  case  of  liquid  sodium  for  future  study.  The  literature  was  searched 
for  correlations  that  would  cover  the  low-to-high  Prandtl  number  range 
as  well  as  the  laminar  and  turbulent  regimes  of  free  and  forced 
convection.  The  correlations  that  were  used  are  presented  below  along 
with  further  references. 


4.2.1  Nusselt  Number  Correlations  for  Free  Convection 

Host  of  the  work  done  in  free  convection  correlations  deals  with 
either  a  constant  temperature  or  constant  heat  flux  situation.  Both 
laminar  and  turbulent  flows  have  been  examined.  Under  laminar  flow 
conditions  Kays  and  Crawford  [57]  state  the  constant  heat  flux  case 
yields  Nusselt  numbers  about  15%  higher  than  the  constant  temperature 
case  while  in  turbulent  flow  both  cases  yield  similar  Nusselt  numbers. 
Since  the  present  work  is  concerned  with  examining  minimum  cooling 
requirements  for  upper  bounds,  the  constant  temperature  correlations 
were  chosen.  Therefore,  if  cooling  conditions  could  be  established  that 
prevented  melting  under  constant  temperature  conditions,  constant  heat 
flux  conditions  would  be  satisfied  too.  Of  course,  the  present  work's 
situation  is  more  like  that  of  an  "unheated  starting  length"  problem, 
but  considering  the  complete  plate  to  be  at  constant  temperature 
provides  a  further  upper  bound  to  the  problem. 

For  readers  interested  in  the  constant  heat  flux  case,  Fuji!  and 
Fuji!  [58]  and  Churchill  and  Ozoe  [59]  examine  the  laminar  flow  regime. 
Turbulent  flow  free  convection  is  examined  by  Vliet  and  Ross  [60]  and 
Vliet  and  Liu  [61].  Though  the  above  listings  are  obviously  not 


comprehensive,  they  will  give  the  interested  reader  an  idea  of  what  work 
has  been  done  in  the  field. 


Nusselt  number  correlations  for  free  convection  are  usually 
functions  of  the  Frandtl  number  and  Grashof  number.  The  Prandtl  number 
of  a  fluid  provides  a  measure  of  the  relative  effectiveness  of  momentum 
transfer  by  diffusion  compared  to  energy  transfer  by  diffusion.  The 
Prandtl  number  is  fixed  by  properties  of  the  fluid  (v:  kinematic 
viscosity  and  a:  thermal  diffusivity)  and  thus  provides  no  latitude  in 
adjusting  the  Nusselt  number.  The  Grashof  number  provides  a  measure  of 
the  ratio  of  buoyancy  forces  to  viscous  forces  in  the  velocity  boundary 
layer.  It  too  is  largely  fixed  by  properties  of  the  fluid  (0: 
volumetric  thermal  expansion  coefficient  and  v)  as  well  as  temperature 
conditions  of  the  problem.  Thus,  as  expected,  when  dealing  with  a  free 
convection  situation  the  main  way  to  adjust  the  Nusselt  number  is  to 
change  the  cooling  fluid.  So  in  examining  cooling  problems  involving 
free  convection,  the  main  question  is  whether  the  cooling  requirements 
can  be  met  with  the  fluids  available.  The  correlations  found  to  answer 
that  question  immediately  follow. 

4. 2. 1.1  Free  Convection* -Laminar  Flow- -Constant  Temperature 

Since  the  present  work  uses  liquid  sodium,  air,  and  water  as  the 
potential  coolants,  correlations  were  needed  for  Prandtl  numbers  in  the 
range  from  approximately  0.005  to  12.0.  For  the  case  of  the  lower  end 
of  this  range,  LeFevre  [62]  determined  limiting  correlations  as  Prandtl 
number  approaches  0  and  infinity.  For  the  present  work,  his  correlation 
for  the  lower  limit  was  used  to  examine  the  case  of  liquid  sodium: 

Nu(x)  -  0. 600(Gr(x)Pr2)0' 25 
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For  the  case  of  air  (Pr  -  0.7),  the  correlation  developed  by 
Churchill  and  Usagi  [63]  was  used: 

Nu(x)  -  0.503Gr(x)1/4Pr1/4/[l  +  (0.492/Pr)9/16]4/9 
For  the  case  of  water  (Pr  -  1.0  and  12.0),  the  expression  reported 
by  Ede  [64]  was  used: 

Nu(x)  -  0. 75[2Pr/5(l  +  2Pr1/2  +  2Pr)]1/4  [Gr(x)Pr]1/4 


4. 2. 1.2  Free  Convection- -Turbulent  Flow-Constant  Temperature 

As  with  laminar  flow,  correlations  applicable  over  the  same  Prandtl 
number  range  were  needed.  For  the  cases  of  liquid  sodium  (Pr  -  0.005) 
and  air  (Pr  -  0.7),  the  work  of  Bailey  [65]  was  employed.  For  sodium, 
the  correlation  used  was: 

Nu(x)  -  0.060Gr(x)1/4 

The  above  correlation  applies  for  the  Grashof  number  range: 

1010  -  1015. 

For  air,  Bailey  found  the  following  applied: 

Nu(x)  -  0.183(Gr(x)Pr)0,31 

This  relation  is  valid  for  air  in  the  Grashof  number  range: 

109  -  1015. 

For  the  case  of  water  (Pr  -  1.0),  the  work  of  Eckert  and  Jackson 

[66]  was  used.  They  found  the  following  applied  for  Pr  -  1.0  and 

Gr(x)  -  1010  -  1012  : 

Nu(x)  -  0.0295[Pr7/(l  +  0.49Pr2/3)6]1/15Gr(x)2/5 
The  reason  the  ranges  of  Grashof  number  in  the  turbulent  flow  regime 
usually  have  a  lower  limit  of  about  Gr(x)  -  10*^  is  because  it  is  around 
this  value  that  turbulent  free  convection  is  realized.  Pitts  and  Sissom 

[67]  state  that  the  transition  from  laminar  to  turbulent  flow  in  natural 


(free)  convection  usually  occurs  when  the  product  of  Prandtl  number  and 

9 

Grashof  number  equals  approximately  10  . 

For  other  work  in  the  area  of  free  convection,  see  the  references 
at  the  end  of  this  work  [68-71]. 


4.2.2  Nusselt  Number  Correlations  for  Forced  Convection 

As  with  free  convection,  much  of  the  work  accomplished  dealing  with 
forced  convection  over  flat  plates  is  concerned  with  either  the  constant 
temperature  case  or  the  constant  heat  flux  case.  Since  Kays  and 
Crawford  [57]  state  that  constant  heat  flux  correlations  yield  Nusselt 
numbers  approximately  36%  higher  in  laminar  flow  and  about  4%  higher  in 
turbulent  flow  than  their  constant  temperature  counterparts,  the 
constant  temperature  correlations  will  again  be  the  ones  used.  (For  the 
constant  heat  flux  case,  Kays  and  Crawford  [57]  also  give  a  summary  of 
the  correlations  available  as  well  as  the  work  done  in  the  field) . 

Also,  like  free  convection,  the  "unheated  starting  length"  problem  is 
being  treated  as  a  constant  temperature  condition  over  the  entire  plate 
in  order  to  establish  an  upper  bound  to  the  problem  of  minimum  cooling 
requirement  determination. 

For  the  forced  convection  correlations,  Nusselt  numbers  are 
expressed  as  functions  of  Prandtl  number  and  Reynolds  number.  The 
Reynolds  number  provides  a  ratio  of  inertia  to  viscous  forces.  Unlike 
its  counterpart  in  free  convection  (the  Grashof  number),  the  Reynolds 
number  can  be  used  to  adjust  the  Nusselt  number  without  having  to  change 
cooling  fluids.  This  adjustment  can  be  made  by  increasing  the  free 
stream  velocity  of  the  fluid  passing  over  the  plate.  Thus,  for  the 
present  study,  the  objective  of  using  the  forced  convection  correlations 
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will  be  to  determine  a  velocity  necessary  so  that  the  cooling  fluid  in 
question  will  meet  the  minimum  cooling  requirements.  If  the  velocity 
determined  is  unattainable  by  normal  processes  used  today,  the  cooling 
fluid  will  be  deemed  "disqualified".  A  summary  of  the  correlations  used 
in  the  present  work  follows. 


4. 2. 2.1  Forced. Convection- -Laminar  Flow- -Constant  Temperature 


Kays  and  Crawford  {57]  present  Nusselt  number  correlations  valid 
for  very  low  to  very  high  Frandtl  numbers.  Their  correlations  are  used 
in  the  present  work  for  the  cases  of  air  and  water.  For  moderate  values 
of  Frandtl  number  (Pr  -  0.5  -  10.0),  they  present  the  following 
correlation: 


Nu(x)  -  0.332Pr1/3Re(x)1/2 

(This  correlation  is  also  reported  by  Incxopera  and  DeWitt  {55).) 

For  very  high  Prandtl  numbers  (Pr  >  10.0),  the  above  correlation  is 
slightly  modified  to: 

Nu(x)  -  0.339Pr1/3Re(x)1/2 

For  the  very  low  Prandtl  number  case,  Eckert  and  Drake  [72]  give  a 
correlation  valid  for  the  range:  0.005  <  Pr  <  0.05.  The  correlation 
used  was: 

Nu(x)  -  [Re(x)Pr]1/2/[1.55Pr1/2  +  3.09(0.372  -  0.15Pr)1/2  ] 

Thus,  this  correlation  was  used  for  examining  liquid  sodium  as  a 
coolant. 

It  is  important  to  note  that  laminar  flow  over  a  flat  plate  has  a 
limiting  Reynolds  number  of  approximately  60,000.  After  this  critical 
Reynolds  number  is  reached,  the  flow  begins  to  transition  to  the 
turbulent  regime. 


i 


For  the  moderate  Prandtl  number  range  (0.5  <  Pr  <  5.0),  valid  for 


air  and  for  water  under  high  pressure,  Kays  and  Crawford  [57]  give  a 
correlation  good  for  Reynolds  numbers  up  to  several  million: 

Nu(x)  -  0.0287PrRe(x)°-8/[0.169Re(x)'°-1(13.2Pr  -  10.16)  +  0.9] 

For  the  case  of  water  under  very  high  pressure  (Pr  -  12.0)  and  the 
case  of  liquid  sodium  (Pr  -  0.005),  Eckert  and  Drake  [73]  give  the 
following  correlation  valid  for  Reynolds  numbers  up  to  10^  : 

Nu(x)  -  0.0296PrRe(x)°‘ 8/[l  +  0.87A(Re(x))‘01(Pr-l)] 


where : 


A 


1.5Pr 


-1/6 


4.3  Cooling  Condition  Determination  and  Coolant  Evaluation  Results 
Vith  the  above  Nusselt  number  correlations,  minimum  cooling 
conditions  can  be  established  and  coolant  suitability  can  be  evaluated. 
In  evaluating  the  minimum  cooling  conditions,  a  minimum  Nusselt  number 
must  be  established  for  each  individual  case.  Thus,  if  cooling 
conditions  can  be  designed  that  exceed  this  minimum  Nusselt  number, 
melting  of  the  metal  plates  will  be  prevented.  To  establish  this 
minimum  local  Nusselt  number,  the  definition  of  Nusselt  number  is 
invoked.  Namely: 

Nu(x)  -  h(x)x/k^ 

The  h(x)  value  is  supplied  by  the  model  and  is  the  heat  transfer 
coefficient  that  will  prevent  any  melting  of  the  plate.  The  distance  x 
in  the  present  work  is  0.025  m  (the  Center  of  the  plate),  since  this  is 
the  location  of  the  highest  temperature  on  the  back  face  (as  well  as  the 


45 


front  face)  of  the  plate. 

The  conductivity  of  the  fluid  is  k^  and  is  evaluated  at  the  average 
temperature  between  the  backside  of  the  plate  and  the  free  stream 
temperature.  The  backside  temperature  of  the  plate  is  supplied  by  the 
model  while  the  free  stream  temperature  is  defined  to  be  300*K.  The 
resulting  minimum  local  Nusselt  numbers  are  given  in  Table  2.  (Note: 

The  special  cases  where  water  under  high  pressure  is  the  coolant  were 
included  for  the  express  purpose  of  examining  a  fluid  with  Prandtl 
numbers  of  1 . 0  and  12.0.) 

Once  the  minimum  cooling  requirements  were  established  for  each 
metal/coolant  c  mbination,  the  Nusselt  number  correlations  could  be 
graphed  to  determine  the  suitability  of  the  coolants  under  forced  and 
free  convection  conditions.  Figures  7-12  show  the  results  for  forced 
convection  correlations.  On  these  graphs,  intersection  of  the  minimum 
local  Nusselt  number  lines  with  the  correlation  lines  yield  minimum 
Reynolds  numbers  necessary  to  prevent  melting.  For  free  convection 
correlations,  Figures  13  and  14  show  the  results.  On  these  graphs, 
intersection  of  the  minimum  local  Nusselt  number  lines  with  the 
correlation  lines  yield  minimum  Grashof  numbers  necessary  to  prevent 
melting. 

4.3.1  Forced  Convection  Results 

One  fact  is  readily  apparent  upon  examining  Figures  7-12;  namely, 
the  laminar  flow  regime  is  totally  disqualified  as  a  possible  cooling 
condition.  This  is  not  really  a  surprise  since  the  problem  does  involve 
a  high  intensity  heat  source  and  since  the  heat  transfer  coefficients 
calculated  by  the  model  were  definitely  in  the  upper  range  of  typical 
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TABLE  2 


Minimum  Local  Nusselt  Numbers 

7  2 

[  0.05  m  Square  Plate,  0.0004  m  thickness,  2  x  10  W/m  Heat  Source  ] 


Metal/Coolant 

Prandtl 

Number 

Average 
Temperature 
(Tfi  +  300)/2 

Conductivity 
of  Coolant 
(kf  :  W/m*K) 

Minimum 
Nusselt  # 
[h(x)0.025/kf] 

Aluminum/ 

Sodium 

0.005 

600°K 

76.40 

9.326 

Stainless  Steel/ 
Sodium 

0.005 

850* 

64.65 

6.845 

Titanium/ 

Sodium 

0.005 

950* 

60.40 

6.167 

Aluminum/ 

Air 

0.7 

600* 

0.0469 

15191.9 

Stainless  Steel/ 
Air 

0.7 

850* 

0.0596 

7424.5 

Titanium/ 

Air 

0.7 

950* 

0.0643 

5793.2 

Aluminum/ 

h2o* 

1.0 

600* 

0.497 

1433.6 

Titanium**/ 

h2o*** 

12.0 

645* 

0.331 

5113.3 

*  H20  under  pressure  of  1.235  x  107  W/m2 

**  Heat  Source  (for  this  case  only)  -  5  x  107  W/m2 

***  H20  under  pressure  of  2.152  x  107  N/m2 


■U'Mum.M  iiwiww  r*  ■■ r  ■  *  ■i*-**'** 
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values  for  forced  convection  liquid  cooling  (according  to  Incropera  and 
DeVitt  [55]).  To  determine  the  velocities  necessary  for  each  coolant, 
each  graph  must  be  examined  separately. 

To  demonstrate  how  the  following  graphs  were  constructed,  the  steps 
taken  to  produce  Figure  9  (the  sodium  coolant  evaluation)  are  displayed 
in  Figures  7*8.  In  Figure  7,  the  minimum  local  Nusselt  number  lines  for 
each  of  the  three  metals  under  construction  are  plotted.  In  Figure  8, 
the  laminar  and  turbulent  Nusselt  number  correlations  are  graphed  for 
forced  convection  conditions  involving  fluids  with  Prandtl  numbers 
between  0.005  and  0.05.  Figures  7  and  8  are  then  combined  to  form 
Figure  9. 

In  Figure  9,  the  turbulent  correlation  line  intersects  the  minimum 
Nusselt  number  lines  for  titanium  and  aluminum  at  Reynolds  number  values 
of  30,000  and  68,000  respectively  (titanium  and  aluminum  are  analyzed  to 
establish  a  range  of  applicable  velocities) .  Rearranging  the  defining 
expression  for  local  Reynolds  number  yields  an  expression  for 
calculating  free -stream  velocity: 

Ufs  “  Re(x)u/x 

Since  x  equals  the  local  distance  in  question  (the  center  of  the  plate, 

0.025  m)  and  v  equals  the  kinematic  viscosity  at  the  average 

-7  2  -7 

temperatures  (approximately  4.0  x  10  m  /s  [aluminum]  and  2.35  x  10 
2 

m  /s  [titanium]),  the  appropriate  free  stream  velocity  can  be 

calculated.  The  result  is  a  free  stream  velocity  of  approximately  0.28 

m/s  for  titanium  and  1.09  m/s  for  aluminum.  Since  these  velocities  are 

easily  attainable,  the  conclusion  is  liquid  sodium  would  serve  as  a 

viable  coolant  for  the  present  situation  of  a  0.0004  m  thick  metal  plate 
7  2 

heated  by  a  2  x  10  W/m  heat  source. 
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In  Figure  10,  the  necessary  Reynolds  number  is  almost  two  orders  of 

magnitude  greater  than  those  in  Figure  9.  Since  the  correlation  being 

graphed  is  only  valid  for  Reynolds  numbers  of  up  to  several  million, 

evaluating  the  case  for  titanium  should  reveal  the  feasibility  of  using 

air  as  the  coolant.  The  correlation  line  intersects  the  titanium 

minimum  cooling  line  at  approximately  a  Reynolds  number  of  5,600,000. 

-  6  2 

The  kinematic  viscosity  of  air  at  950*K  is  112.2  x  10  m  /s.  Again 

using  x  as  the  center  of  the  plate  (0.025m),  the  resulting  free  stream 

velocity  is  approximately  25,133  m/s  (or  56,221  mph) .  Since  this 

velocity  is  unattainable  by  normal  processes ,  air  is  totally 

disqualified  for  forced  convection  cooling. 

In  examining  the  special  case  of  aluminum/water  (Figure  11),  the 

Reynolds  number  at  the  intersection  of  the  minimum  cooling  line  and  the 

Nusselt  number  correlation  is  approximately  780,000.  Using  0.025  m  for 
-7  2 

x  and  1.248  x  10  m  /s  for  the  kinematic  viscosity,  the  free  stream 

velocity  is  approximately  3.89  m/s.  This  result  makes  water  at  1.235  x 
7  2 

10  N/m  a  viable  coolant  for  the  case  of  aluminum. 

The  resulting  Reynolds  number  for  the  special  case  of 

titanium/water  (Figure  12)  is  just  somewhat  below  that  of 

aluminum/water:  Re(x)  -  730,000.  The  conditions  of  the  titanium/water 

problem,  though,  are  completely  different- -with  a  higher  heat  source 

7  2 

intensity  of  5  x  10  W/m  and  a  higher  Prandtl  number  of  12.0.  The 

7  2  -7  2 

kinematic  viscosity  for  water  at  2.152  x  10  N/m  is  1.270  x  10  m  /s. 

Using  0.025  m  for  x  yields  a  free  stream  velocity  of  3.71  m/s.  Thus, 

7  2 

water  at  2.152  x  10  N/m  is  a  viable  cooling  option  for  the  case  of 

7  2 

titanium  subjected  to  a  5  x  10  W/m  heat  source.  Of  course,  the 

7  2 

capability  to  put  the  water  under  a  pressure  of  2.152  x  10  N/m  must 


ALUMINUM 


The  special  case  of  water  under  high  pressure  (1.235  x  10T  N/m2  and 
Pr  -  1.0)  cooling  aluminum  subjected  to  a  2  x  10T  W/m2  localized  heat 
source.  The  minimum  Nusselt  line  for  aluminum  intersects  the 
turbulent  Nusselt  number  correlation  at  a  Reynolds  number  of  about 


TITANIUM 


The  special  case  of  water  under  very  high  pressure  <2.152  x  10T  N/ra2 
and  Pr  -  12.0)  cooling  titanium  subjected  to  a  5  x  10T  W/m2  localized 
heat  source.  The  minimum  Nusselt  number  line  for  titanium  intersects 
the  turbulent  Nusselt  number  correlation  at  a  Reynolds  number  of 
about  730,000. 
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exist  or  water's  viability  in  this  case  is  lost. 

4.3.2  Free  Convection  Results 

Since  air  was  disqualified  as  a  viable  coolant  under  turbulent 

forced  convection  conditions,  it  follows  that  air  is  equally 

disqualified  as  the  coolant  under  free  convection  conditions  since  these 

conditions  have  less  cooling  capability.  Therefore,  a  free  convection 

graph  for  air  is  not  included.  To  determine  what  other  conditions  can 

be  ruled  out  before  graphing,  the  maximum  Grashof  numbers  for  each  case 

need  to  be  determined.  Then  the  product  of  the  Prandtl  number  and 

9 

maximum  Grashof  number  can  be  compared  to  10  to  determine  if  the 

turbulent  regime  needs  to  be  considered.  If  the  product  is  well  below 
9 

10  ,  only  the  laminar  free  convection  correlations  need  be  used. 

In  order  to  calculate  the  maximum  Grashof  number,  it  must  be  broken 
down  into  its  constituent  parts.  Namely: 

Gr(x)  -  g0(Tg  -  TJxV 

where:  g  -  9.807  m/s  (gravitational  constant)  ;  p  -  the  volumetric 

expansion  coefficient  (°K  ;  T  «•  temperature  of  the  backside  of  the 

s 

plate  (°K) ;  T  -  free  stream  temperature  of  the  coolant  (°K)  ;  x  - 

location  on  the  plate  (x  -  0.025  m  for  the  present  work);  and  v  -  the 

2 

kinematic  viscosity  (m  /s) .  The  Grashof  number  will  be  at  a  minimum 

for:  maximum  values  of  p,  minimum  values  of  v,  and  maximum  values  of 

T  -  T  . 
s  ® 

In  examining  sodium  first,  it  is  found  that  titanium  yields  the 
largest  temperature  differences  of  950°.  The  maximum  p  is  approximately 
300  x  10  6  (°K  *■)  and  the  minimum  v  about  2.36  x  10 m^/s.  These 

O 

values  yield  an  estimated  maximum  local  Grashof  number  of  7.84  x  10  and 
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a  Prandtl-Grashof  product  of  3.92  x  10^.  Thus  for  the  sodium  case,  only 
laminar  correlations  are  used. 

For  the  aluminum/water  case  (Pr  -  1.0),  the  local  Grashof  number 
was  calculated  since  there  was  only  one  plate  material.  Thus,  the  value 
obtained  would  also  be  the  maximum  Grashof  number.  For  aluminum, 

T  -  T  -  600°;  x  -  0.025  m;  B  -  5.02  x  10"3(“K-1)  and  v  -  1.248  x 

S  co 

10  7  m2/s.  Therefore,  maximum  local  Grashof  number  -  2.96  x  10*^  which 

also  equals  the  Prandtl-Grashof  product.  Thus,  both  turbulent  and 
laminar  conditions  must  be  considered. 

For  the  titanium/water  case  (Pr  -  12.0),  Ts  -  T^  -  912.5°; 

x  -  0.025  m;  0  -  6.38  x  10“2(*K"1);  and  v  -  1.27  x  10'7m2/s.  Therefore 
the  maximum  local  Grashof  number  *  5.53  x  10^.  Since  the  Prandtl  - 

9 

Grashof  product  is  obviously  well  over  10  ,  a  turbulent  correlation 

needs  to  be  considered  too.  Unfortunately,  only  a  laminar  correlation 

was  found  in  the  literature  for  Pr  >  1.0,  so  the  suitability  of  water  at 
7  2 

2.152  x  10  N/m  as  a  coolant  for  the  titanium  case  had  to  be  estimated 
from  the  aluminum/water  case. 

Thus  with  these  values  of  Grashof  number  and  Grashof- Prandtl 
products,  we  first  decide  to  examine  the  turbulent  correlation  for  the 
aluminum/water  case;  for  if  water  is  found  to  be  an  unsuitable  coolant 
in  the  turbulent  regime,  it  will  necessarily  be  disqualified  in  the 
laminar  regime.  Therefore,  in  Figure  13,  by  noting  where  the  maximum 
Grashof  number  is  plotted  on  the  x-axis,  it  is  readily  apparent  that 
water  is  disqualified  in  the  aluminum/water  case  for  turbulent 
conditions  (and  thus  laminar  conditions).  Also,  for  estimation 
purposes,  the  maximum  Grashof  for  the  titanium/water  case  is  plotted  on 


the  x-axis  and  seen  to  fall  well  below  the  value  needed  for  cooling  in 
even  the  aluminum/water  case.  Assuming  a  similarity  between  the  forced 
convection  and  free  convection  relative  responses  (which  for  forced 
convection  the  Reynolds  numbers  were  within  7%  of  each  other  for  the  two 
cases),  it  is  readily  deduced  that  turbulent  conditions  will  not  yield 
the  necessary  cooling  even  for  the  higher  Prandtl  number  case.  For  the 
aluminum/water  case  alone,  an  increase  of  Grashof  number  of  nearly  25 
times  is  required  to  qualify  water  as  a  viable  coolant.  It  is  very 
unlikely  that  the  titanium/water  case  would  be  much  different. 

Having  disqualified  ater  under  high  pressure  as  a  viable  coolant, 
attention  is  turned  to  Figure  14  where  the  best  chance  for  free 
convection  cooling  conditions  lie.  In  order  to  use  Figure  14  to 
evaluate  sodium  as  a  coolant,  maximum  Grashof  numbers  for  each  material 
have  to  be  calculated  (the  estimated  maximum  Grashof  number  only  defined 
the  flow  regime).  Thus,  for  titanium:  0  *  287.6  x  10*6  (V1);  Ts  - 
-  950*;  and  u  -  2.365  x  10*  m  /s.  Therefore,  maximum  Gr(x)  -  7.49  x 
108.  For  stainless  steel:  0  -  296.6  x  10*6  ('K*1);  T  -  T  -  850* ;  v  - 

S  oo 

-7  2  8 

2.661  x  10  m  /s.  Therefore,  maximum  Gr(x)  -  5.46  x  10  .  And  finally 

for  aluminum:  0  *  270.1  x  10*6  (’K*1);  T  -  T  -  600*;  v  -  3.942  x  10'7 
r  '  s  ® 

2  8 

m  /s.  So,  maximum  Gr(x)  -  1.60  x  10  .  With  these  values  calculated  and 
plotted,  Figure  14  can  now  be  interpreted. 

The  conclusion  readily  evident  from  Figure  14  is  that  sodium  is 
disqualified  as  a  viable  coolant  for  the  aluminum  case.  The  titanium 
case,  on  the  other  hand,  makes  sodium  appear  to  be  promising  as  a 
coolant.  For  the  stainless  steel  case,  sodium  probably  is  not  a  viable 
coolant,  since  a  27.3%  increase  in  Grashof  number  would  be  needed  to 
prevent  melting.  Thus,  we  conclude  that  in  the  laminar  free  convection 


ALUMINUM 
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regime,  only  sodium  is  a  viable  coolant  and  only  for  a  titanium  plate, 

7  2 

0.05  m  square  and  0.0004  m  deep,  heated  by  a  2  x  10  W/m  high  intensity 

-5  2 

heat  source  localized  to  an  8.125  x  10  m  area  in  the  center  of  the 
plate. 

This  conclusion  is  somewhat  surprising  since,  judging  from  the 

typical  values  of  heat  transfer  coefficient  given  by  Incropera  and 

2 

DeVitt  [55],  a  value  of  14,900  W/m  *K  is  well  into  the  upper  ranges  of 
typical  forced  convection  values  in  liquids.  However,  liquid  metals  are 
definitely  atypical,  and  the  viability  of  sodium  as  a  coolant  for  the 
titanium  case  in  free  convection  is  definitely  proof  of  that.  Any 
questions  about  the  accuracy  of  LeFevre's  correlation  [62]  are  put  to 
rest  by  Ede  [64]  in  his  comparision  of  LeFevre’s  correlation  with 
computer  solutions  of  the  analytical  problem.  The  final  result  is  that 
the  free  convection  cooling  possibilities  of  the  titanium/sodium 
combination  are  very  good.  (For  a  general  error  analysis  overview,  see 
Appendix  G . ) 


CONCLUSIONS  AND  FURTHER  STUDY 
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5.1  Conclusions 

From  the  investigation  of  the  basic  problem  of  a  localized  high 
intensity  heat  source  directed  against  one  surface  of  a  thin,  vertical 
metal  plate  with  forced  or  free  convection  cooling  on  the  opposite 
surface,  several  conclusions  can  be  drawn: 


1. 


For  free  convection  cooling,  only  liquid  sodium  is  a  viable 
coolant  and  only  in  the  case  of  a  titanium  plate  subjected  to  a 
2  x  107  W/m2  heat  source. 


Sodium  is  a  viable  coolant  in  the  turbulent  forced  convection 
regime  for  aluminum,  titanium,  and  stainless  steel  subjected  to 
a  heat  source  of  intensity  2  X  107  W/m2.  Velocities  of  no  more 
than  approximately  1  m/s  are  required. 


Water  under  high  pressure  is  also  a  viable  coolant  for  the 
special  cases  of  aluminum  under  an  intensity  of  2  x  107  W/m2 
and  titanium  under  5  x  107W/m2 . 


4. 


Air  is  not  a  viable  coolant  for  either  free  or  forced 
convection  regimes  for  any  of  the  metals  or  intensities 
considered. 


5. 


Heat  transfer  coefficient  as  a  function  of  affected  area  has 
similar  response  curves  for  various  intensities.  In  all  cases, 
there  is  a  rapid  rise  in  the  value  of  the  heat  transer 
coefficient  needed  to  prevent  melting  with  initial  increase  in 
area  and  a  just  as  rapid  leveling  off  after  a  certain 
"critical"  area  is  reached.  This  "critical"  area  occurs  at 
smaller  values  of  affected  area  for  lower  conductivity  metals. 


The  heat  transfer  coefficient  as  a  function  of  heat  intensity 
produces  a  near- linear  response  in  metals  of  high  conductivity 
and  a  near  exponential  response  for  metals  of  low  conductivity. 


The  present  study  appears  to  be  unique  since  the  literature 
revealed  no  similar  studies  in  localized  heat  sources. 


8. 


The  finite  element  program  ANSYS,  by  Swanson  Analysis  Systems, 
Inc.,  is  a  very  accurate  and  versatile  program  that  is 
especially  suited  for  this  study  of  the  basic  problem. 
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Based  on  the  findings  of  the  present  work,  recommended  areas  of 
further  study  are: 


1.  Actual  experiments  to  test  the  findings  of  the  model. 

2.  Investigation  of  the  transient  response  to  a  pulsed  high 
intensity  heat  source  and  the  effect  high  heat  transfer 
coefficients  have  on  such  a  response. 

3.  Varying  the  location  of  the  affected  area  (instead  of  keeping 
it  in  the  center  of  the  plate)  and  noting  the  effect  on  the 
local  minimum  heat  transfer  coefficient. 

4.  Investigation  of  the  effect  that  a  heat  source  with  varying 
intensity  has  on  the  basic  problem. 

5.  Investigation  of  a  moving  or  scanning  heat  source. 

6.  Investigation  of  boiling  as  a  cooling  condition,  and  of  the 
effect  of  certain  surfaces  to  enhance  boiling. 

7.  Examination  of  more  materials  that  span  a  wide  range  of 
conductivities  in  order  to  establish  any  trends  in  the  graph  of 
minimum  heat  transfer  coefficient  versus  conductivity. 

8.  Examination  of  the  effect  on  minimum  heat  transfer  coefficient 
when  radiation  is  taken  into  account. 
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APPENDIX  A 


Sample  Computer  Program  Used  to  Test  Accuracy  of  ANSYS 


C--The  general  analysis  data  generator  is  chosen. 

/PREP7 

/TITLE, TITANIUM  SLAB- -TORVIK  3RD  EX. 

C--The  thermal  analysis  module  is  selected. 

KAN, -1 

C--The  3-D  isoparametric  thermal  solid  element  (STIF  70)  is  chosen. 

C  This  linear-approximation  element  has  8  nodal  points  with 
C  temperature  being  the  single  degree  of  freedom.  The  element  is 
C  applicable  to  a  3-D,  steady-state  or  transient,  thermal  analysis 
C  [52]. 

ET, 1 , 70 

C- -Conductivity  entered. 

KXX.1,0.145 
C- -Density  entered. 

DENS, 1,4. 43 

C- -Specific  heat  entered. 

C, 1,0. 77 

C--Node  pattern  is  generated 
N ,  1 

N.21,2.5 

FILL 

NGEN, 21, 21, 1,21,1, ,0.125 
NGEN, 5, 500, 1,450,1, , ,0.01 
C- -First  element  is  defined. 

E, 1,2, 23, 22, 501, 502, 523, 522 

C--Element  sets  are  generated  using  the  previously  defined  element(s). 
EGEN , 20 , 1 , - 1 
EGEN.4 , 500 , -2C 
EGEN, 20, 21, -80 

C--The  absolute  value  of  the  first  number  after  ’ITER,'  is  the  maximum 
C  iterations  for  the  first  load  step.  ’TIME'  -  time  at  end  of  first 
C  load  step.  Thus,  the  first  interval  -  0.0964  -  0.0  and  the  intera- 
C  tion  time  step  (ITS)  -  TIME/  1-160  I.* 

ITER, -160,0, 160 
TIME, 0.0964 

C- -Loads  are  set  up  as  step  loads. 

KBC.l 

C- -Initial  temperature  is  defined. 

TUNIF, 300 

C- -Analysis  is  set  up  to  terminate  if  solution  does  not  converge. 
CNVR.0.1, , , ,1 

C- -Approximate  Gaussian  distribution  applied  to  an  approximate 

circular  area  of  0.005m  radius.  The  distribution  is  of  the  form: 

F(r)  -  F  x  exp(-r^/2o)  ;  where  a  -  0.0025  m 


The  values  in  the  program  were  calculated  by  taking  the  radial 
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distance  of  the  center  of  the  face  of  each  affected  element  and 
substituting  that  distance  into  the  equation.  The  resulting 
values  are  the  approximate  average  intensities  for  the  affected 
elements .  The  exception  is  the  central  element  which  uses  the  peak 
intensity  due  to  the  nature  of  the  Gaussian  distribution. 
CV,1,2,2.0E-5,1E8,  ,  ,23,22 
CV, 22 , 23 , 1 .462E-5 , 1E8 , , ,44,43 
CV,43,44,8.86E-6,1E8, , ,65,64 
CV,64,65,4.19E-6,1E8, , ,86,85 
CV,2,3,1.462E-5,1E8, , ,24,23 
CV,3,4,8.86E-6,1E8, , ,25,24 
CV,4,5,4.19E-6,1E8, , ,26,25 
CV,23,24,1.14E-5,1E8, , ,45,44 
CV,24,25,6.93E-6,1E8, , ,46,45 
CV, 25 , 26 , 3 . 26E-6 , 1E8 , , ,47,46 
CV,44 ,45 , 6 . 93E-6 , 1E8 ,,,66,65 
CV,45,46,4.19E-6,1E8, , ,67,66 
CV, 65 , 66 , 3 . 26E-6 , 1E8 ,,,87,86 

C--The  nodes  are  listed  from  which  temperature  printouts  are  required. 
PRDISP.1,1,5,1 
,1.22,26,1 
.1,43,47,1 
,1,64,67,1 
,1,501,2001,500 

C- -Stress  data  for  the  last  iteration  is  printed. 

PRSTR, 160 

C- -Stress  data  for  the  last  iteration  is  saved  for  post-processing  and 
C  written  to  File23. 

POSTR, 160 
LWRITE 

C--Second  load  step  is  defined.  Interval  -  .1042-. 0964 
C  ITS  -  13/( . 1042- .0964) 

ITER, 13, 1,1 
TIME, 0.1042 

C--Data  written  to  File23. 

LWRITE 

C- -Third  load  step  defined.  Interval  -  .1064-. 1042 
C  ITS  -  1-2  l/(.1064-.  1042) 

ITER, -2,1,1 
TIME, 0.1064 

C--Data  written  to  File23. 

LWRITE 

C- -Fourth  load  step  defined.  Interval  -  .1064-. 1042 
C  ITS  -  1-6  l/(.ll  -  .1064) 

ITER, -6,1,1 
TIME, 0.11 

C- -Stress  data  for  last  iteration  printed  and  saved 
PRSTR, 6 
POSTR, 6 

C--Data  written  to  File23 
LWRITE 

C- -Complete  data  files  written  to  File27. 

AFWRIT 

C-Exit  general  analysis  data  generator. 
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FINISH 

/EXEC 

C- -Commence  solution  phase  with  File27  data. 

/INPUT, 27 

FINISH 

/POST26 

C--Record  all  data  for  the  following  nodes:  1-5,  22-26,  43-47,  64-68, 
C  85-89,  501,  1001,  1501,  2001. 

NUMVAR, 30 
DISP, 2,1, TEMP 
DISP, 3 , 2 , TEMP 
DISP, 4, 3, TEMP 
DISP, 5, 4, TEMP 
DISP, 6, 5, TEMP 
DISP, 7, 22, TEMP 
DISP, 8, 23, TEMP 
DISP, 9, 24, TEMP 
DISP, 10, 25, TEMP 
DISP, 11, 26, TEMP 
DISP, 12, 43, TEMP 
DISP, 13, 44, TEMP 
DISP, 14, 45, TEMP 
DISP, 15, 46, TEMP 
DISP, 16, 47, TEMP 
DISP, 17, 64, TEMP 
DISP, 18, 65, TEMP 
DISP, 19, 66, TEMP 
DISP, 20, 67, TEMP 
DISP, 21, 68, TEMP 
DISP, 22, 85, TEMP 
DISP, 23, 86, TEMP 
DISP, 24, 87, TEMP 
DISP, 25, 88, TEMP 
DISP, 26, 89, TEMP 
DISP, 27, 501, TEMP 
DISP, 28, 1001, TEMP 
DISP, 29, 1501, TEMP 
DISP, 30, 2001, TEMP 

C- -Printout  the  data  stored  over  the  time  interval  0.06  -  0.11 
C  seconds  in  the  following  groupings  of  nodes  (node  numbers  follow 
C  PRVAR) . 

PRTIME, 0.06, 0.11 
PRVAR, 2, 3, 4, 5, 6 
PRVAR, 7, 8, 9, 10, 11 
PRVAR, 12, 13, 14, 15, 16 
PRVAR, 17, 18, 19, 20, 21 
PRVAR, 22, 23, 24, 25, 26 
PRVAR, 2 ,27,28,29,30 
FINISH 
/EOF 


♦Note:  In  transient  thermal  analysis  for  ANSYS,  the  integration  time  step  (ITS) 

is  related  to  the  "conducting  length"  (S)  of  an  element  in  the  region 
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over  which  the  largest  gradient  acts  (usually  the  length  parallel  to  the 
heat  flow).  The  larger  the  thermal  gradient,  the  smaller  both  the  ITS 
and  6  should  be.  The  following  guideline  is  given  to  help  the  user 
select  an  initial  ITS  for  a  given  model:  (ITS).  £  S2/Ua  where  (ITS),  is 
the  initial  ITS  and  a  is  the  material  thermal  aiffusivity.  Smaller1ITS 
sizes  (or  larger  elements)  may  cause  temperature  oscillations  in  the 
large  temperature  gradient  region. 


APPENDIX  B 


Sample  Computer  Program  for  Determining  Minimum 
Heat  Transfer  Coefficient 

C--The  general  analysis  data  generator  is  chosen. 

/PREP7 

C--New  heat  transfer  coefficient  is  selected. 

*SET, FILM, 1.49 

C- -Title  is  changed  to  reflect  new  values. 

/TITLE,  TITANIUM  SLAB  --  FILM  -  1.49 
C- -Thermal  analysis  module  is  selected. 

KAN, -1 

C- -Choice  of  element  type  (STIF  70)  is  made. 

ET,1,70 

C- -Temperature  table  of  thermal  properties  is  entered.  Since  more  than 
C  6  values  a  2nd  line  is  needed. 

MPTEMP , 1 , 100 , 200 , 400 , 600 , 8 00 , 1000 
MPTEMP , 7 , 1200 , 1500 

MPDATA.KXX, 1,1, 0.305, 0.245, 0.204, 0.194, 0.197, 0.207 

MPDATA.KXX, 1 , 7 , 0. 220, 0. 245 

C--Node  pattern  is  now  established. 

N.l 

N,21,2.5 

FILL 

NGEN.21, 21, 1,21,1, ,0.125 
NGEN, 5, 500, 1,450,1, , ,0.01 
C- -First  element  is  defined. 

E, 1,2, 23, 22, 501, 502, 523,522 

C--More  sets  of  elements  are  generated  using  the  previous  set(s)  as  a 
C  pattern. 

EGEN.20,1, -1 
EGEN.4, 500, -20 
EGEN, 20, 21, -80 

C- -Maximum  of  ten  iterations  is  set. 

ITER, -10,0,10 

C--Load  is  defined  as  a  step  load. 

KBC.l 

C- -Initial  temperature  for  material  properties  is  defined. 

TUNIF, 300 

C- -Convergence  criteria  of  0.1  degrees  is  set. 

CNVR.0.1 

C--Back  surface  of  plate  is  assigned  the  heat  transfer  coefficient 
C  defined  at  the  start  of  the  program. 

CVSF, ,3, 0.04, FILM, 300  ?  2 

C- -Front  surface  of  plate  has  a  2  x  10  W/ni  intensity  applied  to  a 
C  limited  area. 

CV , 1 , 2 , 2E-5 , 1E8 ,4 , 1,23,22 
CV,22,23,2E-5, 1E8 ,25,1,44,43 
CV,43 ,44 , 2E-5, 1E8 ,45 ,1,65,64 
CV, 64 , 65 ,2E-5 , 1E8 ,65,1,86,85 
C- -Analysis  is  written  to  File27. 

AFWRIT 
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C- -General  analysis  data  generator  is  exited. 

FINISH 

C- -Solution  phase  is  started  with  File27  data. 

/EXEC 
/INPUT, 27 
FINISH 

C--The  ten  highest  temperatures  are  printed  out  along  with  their 
C  coordinates. 

/P0ST1 
SET,  1 

NSORT , TEMP , , ,10 
PRTEMP 
PRNODE 
FINISH 

C--ANSYS  is  exited. 

/EOF 
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APPENDIX  C 

Minimum  Heat  Transfer  Coefficients  Necessary  to  Prevent  Melting  at 

Various  Intensities 
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1  x  107 

10,500 

6400 

7600 
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28,500 

14,900 

17,700 

3  x  107 

47,800 

25,700 

30,400 

4  x  107 

68,100 

40,600 

48,000 

5  x  107 

89,700 

67,700 

81,500 
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130,200 

165,200 
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APPENDIX  D 
DATA  FOR  FIGURES  2-4 

(Note:  In  the  listings  below,  the  data  for  the  lines  graphed  in  Figures 

2-4  are  presented.  If  the  data  point  was  calculated  exactly  by  the 

model,  the  number  stands  alone.  Otherwise,  [  ]  indicates  the  value  was 

extrapolated  from  model  generated  data;  (  )  indicates  the  value  was 

interpolated  from  model  generated  data;  (  )  indicates  the  value  was 

estimated  from  studying  curve  behavior.  The  appropriate  data  for  each 

intensity  is  listed  underneath  the  respective  intensity  across  from  the 

appropriate  affected  area.  The  left  column  is  the  area  affected  by  the 
2  -4 

heat  source  in  m  x  10  .  The  values  on  the  right  underneath  the 

2  4 

intensities  are  heat  transfer  coefficients  in  W/m  °K  x  10  .  Though  not 
listed,  the  0.0  heat  transfer  coefficient  for  0.0  area  is  a  point  for 


each  line.) 


Figure  2: 

ALUMINUM 

Area 

1  x  107  W/m2 

2  x  107  W/m2 

1 

x  107  W/m 

0.25  : 

(0.26) 

1.46 

(3.30) 

0.5  : 

(0.72) 

2.40 

{4.46} 

0.8125  : 

1.05 

2.85 

4.78 

1.625  : 

(1.41} 

3.20 

(5.04) 

2.4375  : 

(1.53} 

3.30 

{5.14} 

3.25  : 

(1.57) 

3.33 

(5.18) 

4.0625  : 

1.59 

3.34 

(5.20) 

4  x  107  W/m2 

5  x  10 7  W/m2 

L 

x  107  W/m1 

(5.20) 

(7.23) 

(9.34) 

(6.45) 

(8.57) 

[10.83] 
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1.625  : 

{0.78} 

1.79 

(3.06) 

2.4375  : 

{0.79} 

1.80 

(3.07) 

3.25  : 

0.80 

1.80 

{3.08} 

4  x  IQ7  W/m2 

5  x  10 7  W/m2 

6  X  107  U/m 

{4.57} 

{7.46} 

(14.47) 

{4.76} 

(8.03) 

{16.16} 

4.80 

8.15 

16.52 

(4.84) 

{8.28} 

{16.98} 

{4.89} 

{8.41} 

(17.44} 

{4.93} 

{8.54} 

[17.9] 
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APPENDIX  E 

Minimum  Heat  Transfer  Coefficient  as  a  Function  of  Thickness 

7  2 

[Results  from  a  2  x  10  W/m  Local  High  Intensity  Heat  Source] 
Aluminum  Thickness  Heat  Transfer  Coefficient 


(m) 

(W/m2  °K) 

(0.0) 

(31,600)* 

0.0001 

31,700 

0.0002 

30,900 

0.0004 

28,500 

0.0008 

23,100 

0.0016 

14,500 

0.0032 

6,200 

*  Note:  31,600  is  the  theoretical  heat  transfer  coefficient  for  a 
plate  of  'zero'  thickness  (i.e.  100%  heat  conduction). 
[Included  for  comparison.] 
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APPENDIX  F 


Data  for  Heat  Transfer  Coefficients  Versus  Conductivity 

7  2 

[Values  listed  are  for  a  heat  source  of  2  x  10  W/m  ] 


Conductivity 

Heat  Transfer  Coefficient 

(W/m  ®K) 

(W/m2  *K) 

24.5* 

Titanium 

14,900 

31.7 

Stainless  Steel 

17,700 

218 

Aluminum 

28,500 

*  Note:  Values  of  conductivity  are  those  used  by  the  mode-1  taken  from 
the  table  of  values  from  Incropera  and  DeWitt  [55].  In  each 
case,  the  conductivity  was  based  on  the  average  temperature 
through  the  plate- -which  resulted  in  using  the  highest  value 
in  each  table. 
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APPENDIX  G 

Error  Analysis  Overview 

In  examining  the  sources  of  error  in  this  study,  the  two  major 

sources  can  be  identified  as:  the  computer  model  itself,  and  the  errors 

in  the  correlations  used  to  evaluate  the  coolants .  As  noted  in  Table  1 , 

the  model  yields  an  initiation  of  melting  time  only  1.45%  below  the 

exact  analytical  solution  (computer  generated) .  Since  the  model  uses 

square  elements  to  approximate  a  circular  area  and  since  the  area  of  the 

model  is  3.45%  greater  than  that  of  the  analytical  solution,  this  lower 

melting  time  is  expected.  The  conclusion  then  is  that  the  error  the 

model  incurs  when  examining  the  metal  plates  heated  by  localized  heat 

7  2 

sources  in  the  range  of  1  -  6  x  10  W/m  appears  to  be  approximately  1%. 

One  cautionary  note  is  in  order  for  using  the  model  to  analyze 
small  areas  at  high  intensities.  Namely,  as  the  affected  area 
decreases,  the  number  of  elements  available  for  resolving  the 
temperature  distribution  decreases  proportionately.  Thus,  if  this 
decreased  area  condition  is  coupled  with  an  increase  in  intensity, 
errors  wil^  be  induced  at  the  higher  gradients  since  using  too  few 
elements  to  analyze  a  high  temperature  gradient  situation  usually  yields 
excessively  high  temperatures.  This  is  because  of  the  linear 
approximations  used  in  the  finite  element  program.  Therefore,  if  faced 
with  a  small  area/high  intensity  problem,  the  model  must  have  the  number 
of  elements  covering  the  plate  increased,  at  least  in  the  vicinity  of 
the  area  of  the  plate  subjected  to  the  localized  heat  source.  For  the 
regular  version  of  ANSYS,  this  is  no  problem.  For  the  educational 
version,  with  its  limited  memory  and  wave  front,  increasing  the  elements 


for  this  problem  is  out  of  the  question  unless  the  portion  of  the  plate 

being  examined  is  decreased  (i.e.  invoke  symmetry  to  reduce  from  1/4  to 

1/8  of  the  plate).  This,  indeed,  became  necessary  to  obtain  data  for 

7  2 

intensities  in  the  range  of  5  -  6  x  10  W/m  coupled  with  small  affected 
areas  for  the  cases  of  titanium  and  stainless  steel.  Aluminum's  high 
conductivity  made  this  procedure  unnecessary. 

In  examining  the  errors  in  the  correlations  used  to  evaluate  the 
coolants,  the  literature  indicates  most  to  be  under  10%  (well  under  in 
many  of  the  laminar  cases).  Those  approaching  20%  disagreement  (with 
some  of  their  applicable  experimental  data)  were  for  air,  and  air  was 
disqualified  as  a  viable  coolant  under  both  forced  and  free  convection 
conditions.  Even  if  we  remain  conservative  and  assume  a  15%  error 
possible  in  the  correlations,  the  combination  of  errors  from  the  model 
and  from  the  correlations  still  equals  practically  the  correlation 
error.  Thus,  the  coolant  evaluations  accomplished  in  this  work  should 
not  be  in  error  by  much  more  than  the  uncertainty  of  the  correlation 
used  in  each  particular  Instance.  In  any  case,  for  the  correlations 
used  in  this  study,  no  error  should  be  greater  than  15%.  And,  in  many 
instances,  the  error  will  be  well  under  10%. 

Therefore,  for  any  user  wishing  to  utilize  this  computer  model  to 
determine  minimum  heat  transfer  coefficients  but  supply  his  own 
correlation  to  evaluate  cooling  conditions,  he  can  safely  assume  that 
any  resultant  errors  equal  approximately  the  uncertainty  in  the 
correlation. 
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